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ABSTRACT 

Under  this  research  grant,  the  framework  and  tools  were  developed  for  model  reduction  of 
atomic-scale  many  body  systems.  The  state-affine  mathematical  structure  of  the  original  and 
reduced-order  models  enabled  the  implementation  of  control  through  dynamic  programming  and 
estimation  through  an  extended  Kalman  filter.  In  the  framework  developed  here,  the  state  of  the 
high-dimensional  stochastic  system  is  first  quantified  using  a  high-dimensional  pair  correlation 
function.  This  state  is  then  reduced  using  linear  and  nonlinear  principal  component  analysis  and 
is  discretized  using  self-organizing  maps.  To  create  the  dynamic  model,  a  cell  map  is 
constructed  using  short  simulations  to  quantify  input-dependent  transitions  between  the  discrete 
states.  The  error  associated  with  the  model  reduction  was  quantified  and  analyzed,  and  a  method 
for  predicting  this  error  was  proposed. 

Specific  applications  in  materials  processing  were  considered,  which  motivated  and  guided  the 
development  of  the  model  reduction  framework  and  tools.  A  model  of  gallium  arsenide 
deposition  was  used  to  demonstrate  the  model  reduction  framework.  A  second  modeling  study- 
in  the  molecular  architecture  of  hyperbranched  polymers  was  performed,  and  enabled  a 
comparison  of  common  themes  and  system  specific  features  between  the  two  different 
applications. 
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I.  OBJECTIVES  AND  SIGNIFICANCE 


The  goal  of  this  research  program  was  to  develop  a  systematic  control  methodology  for  state- 
affine  probabilistic  systems,  This  class  of  systems  represents  local  interactions  among  a  large 
number  of  discrete  particles,  and  is  particularly  important  for  dynamics  at  the  nanometer  scale, 
where  the  interactions  among  many  atoms,  molecules,  or  spins  require  a  probabilistic 
description.  Stochastic  realizations  may  be  performed  using  Monte  Carlo  simulations,  but  it  is 
difficult  to  use  them  directly  in  standard  algorithms  for  dynamic  optimization  and  controller 
design.  The  three  main  objectives  of  this  research  program,  as  articulated  in  the  original 
proposal,  were 

1.  to  develop  the  theory  underlying  a  model  reduction  idea  for  kinetic  Monte  Carlo 
simulations 

2.  to  develop  systematic  tools  for  the  reduction  process  and  for  control  of  the  reduced 
models 

3.  to  investigate  the  applicability  of  the  theory'  and  procedures  in  several  application  areas  of 
interest  for  communications,  computation,  and  materials  development. 

Dynamic  models  based  on  physical  principles,  such  as  conservation  of  mass,  momentum,  and 
energy,  are  routinely  used  in  the  design  of  systems  and  controllers.  However,  models  of  atomic 
scale  phenomena,  which  often  are  high  dimensional,  nonlinear,  and  stochastic,  are  not  typically 
used  in  design.  Air  Force  systems  demand  new  high  performance  materials  and  electronic 
devices,  and  due  to  the  myriad  possible  design  and  processing  options  for  nanostructured 
materials,  models  are  needed  to  design,  optimize,  and  control  the  processes.  One  illustrative 
example  is  the  heterostructure  field  effect  transistor  under  study  at  the  Air  Force  Research 
Laboratory  at  Wright  Patterson  Air  Force  Base.  Professor  Gallivan  spent  the  summer  of  2006  at 
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AFRL  working  with  a  team  there  to  model  the  degradation  dynamics  of  these  new  electronic 
devices  for  high-speed  high-power  radar  applications.  Although  these  nanostructured  transistors 
initially  perform  well,  they  degrade  over  time.  It  is  believed  that  the  degradation  and  failure 
modes  depend  strongly  on  the  processing  conditions  and  initial  material  defects.  Due  to  the 
small  scales  of  the  device,  the  models  require  a  consideration  of  a  small  number  of  point  defects 
at  random  locations  in  the  device. 

Two  additional  applications  were  considered  in  this  AFOSR-sponsored  project.  A  lattice  Monte 
Carlo  simulation  of  gallium  arsenide  surface  processing  was  used  to  motivate,  guide,  and 
demonstrate  the  model  reduction  framework  developed  here.  Gallium  arsenide  is  a  ill-V 
compound  semiconductor  used  in  high-speed  transistors,  and  its  interface  properties  must  be 
tightly  controlled  in  layered  electronic  devices  such  as  transistors.  Hyperbranched  polymers  are 
the  second  application  considered  in  this  program.  Their  unique  molecular  architecture  makes 
them  promising  materials  for  sensors  and  catalysis,  although  the  relationship  between  the  process 
inputs  and  the  final  molecular  structure  has  not  previously  been  well  understood. 

The  modeling  framework  proposed  and  developed  here  has  wide  potential  applicability  to  a 
range  of  materials  processing  applications  that  are  modeled  by  stochastic  many  body 
simulations.  The  framework  and  tools  draw  on  a  range  of  approaches  in  the  model  reduction  and 
data  mining  disciplines.  The  model  building  approach  requires  significant  computation,  but  it 
can  be  efficiently  parallelized,  and  the  resulting  reduced-order  model  requires  minimal 
computation,  which  is  critical  for  real-time  controller  implementation.  The  model  building 
process  is  largely  automated  by  the  toots  described  here,  although  some  user  input  is  required  to 
select  an  appropriate  pair  correlation  function  to  define  the  original  state  space.  In  this  report. 
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the  modeling  framework  and  tools  are  first  presented,  followed  by  a  demonstration  of  the 
modeling  and  control  strategy,  after  which  the  applications  are  discussed. 
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Figure  I.  Reduced-order  modeling  framework  and  tools  for  atomic  simulations. 


n.  SIGNIFICANT  ADVANCEMENTS 

The  concept  behind  the  model  reduction  framework  is  illustrated  by  Figure  I.  The  labels  on  the 
arrow's  of  Figure  I  give  examples  of  tools  that  can  be  used  to  implement  each  of  the  steps.  Some 
of  the  concepts  behind  this  reduction  approach  were  previously  proposed  and  used  by  Gal  I  i  van 
and  Murray  [  1 1.  A  goal  of  the  work  supported  by  this  grant  was  to  propose  a  complete  reduction 
process,  and  to  identify  and  demonstrate  tools  that  could  be  used  to  automate  the  reduction. 
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1.  Quantification  of  the  state  space 

A  first  critical  step  in  constructing  a  dynamic  model  is  to  define  the  state  space  of  the  system, 
which  is  a  non-trivial  task  for  an  atomic-scale  many  body  simulation.  Although  the  state  of  the 
stochastic  simulation  can  be  defined  as  the  position  of  each  atom  in  the  sy  stem,  each  realization 
of  the  stochastic  simulation  will  be  expected  to  yield  a  completely  different  state.  However,  the 
control  objectives  for  the  system  are  based  on  overall  statistical  measures  that  define  material 
structure  and  properties,  and  this  is  not  expected  to  change  significantly  between  realizations  that 
are  performed  under  the  same  nominal  process  inputs. 

One  way  to  think  about  this  issue  is  to  consider  symmetries  in  the  stochastic  simulations. 
Periodic  boundary  conditions  are  used,  so  translational  symmetries  exist.  Thus,  it  is  not  the  exact 
position  of  each  atom  that  is  relevant  to  the  state,  but  the  relative  positions  of  the  atoms. 
Furthermore,  most  of  the  atoms  in  the  simulations  are  indistinguishable:  in  the  GaAs 
simulations,  there  are  two  atom  types:  gallium  and  arsenic;  in  the  polymer  simulations,  there  are 
tw'o  types  of  molecules,  which  we  refer  to  as  A:  and  B3.  However,  in  the  entire  simulation 
domain  there  are  thousands  or  millions  of  each  atom.  Thus,  the  state,  especially  in  a  reduced- 
order  model,  need  not  contain  the  exact  position  of  each  atom,  but  instead  should  include  only 
the  relative  positions  of  each  atom  type.  This  concept  is  already  used  in  analyzing  atomic  scale 
simulations,  through  the  pair  correlation  function.  It  is  still  a  high-dimensional  description,  but 
captures  the  frequency  of  pairs  of  atoms  as  a  function  of  their  distance  and/or  angle.  In  surface 
simulations,  such  as  in  the  GaAs  system,  a  height-height  correlation  function  is  often  used  to 
quantify  the  simulation  results,  and  is  a  pair  correlation  function  based  on  the  height  of  the 
surface  at  different  distances  in  x-  and  y-directions  that  define  the  surface.  This  height-height 
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correlation  function  can  in  fact  be  computed  using  the  Fourier  transform  of  the  surface  height,  so 
one  interpretation  of  a  pair  correlation  function  is  in  analogy  to  taking  a  Fourier  transform  of  the 
simulated  surface,  for  the  purpose  of  removing  symmetries.  In  our  GaAs  demonstration,  we 
instead  use  a  pair  correlation  function  based  on  the  distance  between  steps  in  the  surface,  which 
is  instead  more  like  the  Fourier  transform  of  the  derivative  of  the  surface  height.  Due  to  the 
small  changes  in  surface  height  due  to  atomic-height  steps,  we  found  that  the  height-height 
correlation  function  was  extremely  sensitive  to  noise  in  the  simulations,  while  the  step-step 
correlation  function  better  captured  the  discrete  changes  due  to  atomic  height  steps  using  a  small 
number  of  modes.  Thus,  the  selection  of  the  best  pair  correlation  function  for  a  particular  system 
cannot  be  fully  automated  at  this  time,  and  requires  selection  based  on  an  understanding  of  the 
key  material  structures  in  the  system.  Future  work  could  involve  a  more  automated  selection 
process,  although  it  should  also  be  expected  that  some  understanding  of  the  physical  system  will 
aid  in  the  modeling  and  reduction  process. 

3.  Tools  for  model  reduction 

A  significant  accomplishment  of  this  grant  is  the  identification  of  a  set  of  tools  that  can  automate 
the  model  reduction  process  in  a  computationally  tractable  manner.  These  tools  include  principal 
component  analysis  (PCA),  nonlinear  principal  component  analysis  (NLPCA),  self-organizing 
map  (SOM),  cell  mapping,  and  k-nearest  neighbors  interpolation. 

Once  the  pair  correlation  function  has  been  identified  for  the  atomic  simulation,  methods  are 
required  to  reduce  its  dimension  in  order  to  construct  a  reduced-order  model.  The  dimension  of 
the  pair  correlation  function  depends  on  the  simulation  size,  and  should  be  expected  to  scale  with 
either  the  length  or  volume  of  the  simulation  domain,  depending  on  whether  the  correlation 
function  includes  distances  only,  or  also  directions.  However,  beyond  a  certain  correlation 
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length,  the  correlation  functions  should  reach  a  plateau  value  (this  is  a  requirement  in  the  atomic 
simulations  to  ensure  that  the  simulation  domain  is  sufficiently  large).  Additional  redundant 
features  often  exist  in  these  correlation  functions,  which  can  be  associated  with  dominant 
material  structures.  As  a  result,  the  next  step  in  the  modeling  framework  is  to  search  for  lower- 
order  behavior  in  the  molecular  simulations,  using  the  method  of  snapshots  commonly  used  for 
model  reduction  of  partial  differential  equations,  such  as  in  fluid  flow  [2].  In  this  approach, 
simulations  are  run  under  a  wide  range  of  inputs  and  initial  conditions,  and  principal  component 
analysis  (also  called  proper  orthogonal  decomposition)  is  used  to  identify  a  lower  order  state 
space  spanned  by  orthogonal  basis  vectors.  The  coordinates  in  this  new  basis  are  then  used  as 
the  coordinates  in  a  new  reduced-order  model,  Nonlinear  correlations  between  the  states  may 
also  exist,  and  we  have  demonstrated  the  use  of  nonlinear  principal  component  analysis,  which  is 
most  efficiently  performed  after  linear  correlations  have  been  removed  using  principal 
component  analysis  via  autoassodative  neural  nets  (Gal li van  in  2004.  Proceedings  of  DYCOPS). 

A  further  step  we  take  in  our  modeling  method  is  to  discretize  this  low  order  state,  which  enables 
construction  of  a  Markov  chain  model  based  on  input-dependent  transitions.  We  have 
accomplished  this  discretization  using  self-organizing  maps  (SOM)  [3],  although  other  methods 
for  clustering  data  could  alternatively  be  used.  The  self-organizing  map  converges  with  minimal 
computation  required,  and  has  the  additional  beneficial  feature  of  creating  a  projection  of  the 
nonlinear  space  onto  a  two-dimensional  plane  (by  organizing  the  clusters),  which  can  aid  in 
visualization  of  trajectories  of  the  system.  Our  overall  approach  to  state  reduction  was  published 
by  Oguz  and  Gall  i  van  in  2005  in  the  International  Journal  of  Nonlinear  and  Robust  Control. 

Of  course,  it  is  also  possible  to  not  discretize  the  state  space,  and  instead  to  construct  a  low-order 
dynamic  model  on  the  continuous  space  after  performing  the  principal  component  analysis.  This 
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is  accomplished  in  applications  such  as  fluid  flow  using  a  Galerkin  projection  on  the  original 
system  of  partial  differential  equations.  However,  in  the  atomic  simulations,  there  is  no  PDE  for 
the  projection,  and  it  is  not  clear  what  should  be  the  mathematical  form  of  reduced-order 
dynamic  model.  One  option  is  to  project  onto  the  probabilistic  master  equation  that  defines  the 
probability  distribution  of  the  entire  system  over  all  possible  realizations  on  the  stochastic 
simulations.  However,  due  to  the  extremely  large  or  even  infinite  number  of  possible 
configurations  of  the  system,  direct  mathematical  operation  on  this  master  equation  is 
prohibitive. 

For  this  reason,  we  have  chosen  to  discretize  the  reduced-order  state  space  and  construct  a 
dynamic  model  based  on  input-dependent  transitions  between  the  discrete  states  in  the  model. 
As  with  any  reduced-order  model  based  on  the  method  of  snapshots,  it  is  not  expected  that  this 
model  will  be  able  to  predict  states  that  are  significantly  different  from  those  in  the  original  set  of 
snapshots;  however  this  is  not  a  particular  limitation  of  the  discrete  state  space,  but  rather  of  the 
method  of  snapshots.  The  physical  interpretation  of  each  discrete  state  is  that  it  represents  a 
typical  type  of  surface  or  material  configuration  that  was  observed  during  the  method  of 
snapshots  training.  Similar  snapshots  are  grouped  by  the  SOM  into  equivalence  classes,  so  that 
each  snapshot  is  associated  with  one  of  the  discrete  states.  Conversely,  each  discrete  state  is 
associated  with  a  point  in  the  reduced  continuous  state  space,  and  therefore  is  also  associated 
with  the  pair  correlation  function  that  can  be  reconstructed  from  its  point  in  the  reduced 
continuous  state  space. 

Once  the  discrete  state  space  has  been  defined,  the  dynamic  model  is  constructed  by  running 
short  simulations  for  a  particular  predefined  interval  (typically  this  interval  would  be  based  on 
time).  A  simulation,  or  a  set  of  realizations,  is  run  beginning  with  a  configuration  in  each 
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discrete  state,  and  is  run  for  the  discretized  time  step  under  a  particular  process  input.  Here  we 
discretize  the  inputs  into  discrete  levels,  and  construct  transitions  for  each  of  these  levels 
(Interpolation  between  these  levels  in  discussed  in  the  next  section).  At  the  end  of  each 
simulation,  the  pair  correlation  function  and  corresponding  PCA  reduced  coordinates  are 
computed,  and  are  matched  to  the  discretized  state  (or  states)  that  best  matches,  based  on  a 
Euclidean  distance  of  the  PCA  coordinates.  By  performing  short  simulations  beginning  in  each 
discrete  state,  and  at  each  input  level,  a  dynamic  model  is  constructed.  The  computational 
requirements  of  this  reduction  depend  on  the  number  of  discrete  states,  the  number  of  input 
levels,  and  the  number  of  realizations  performed  for  the  purpose  of  averaging  out  random  noise. 
In  our  GaAs  demonstration,  we  use  approximately  150  states,  8  input  levels,  and  10  realizations 
of  each  run.  Since  each  realization  takes  approximately  l  hour,  the  computational  time  on  a 
single  processor  is  significant  but  not  intractable.  However,  the  realizations  can  be  distributed 
across  multiple  processors,  which  thus  reduces  the  time  by  a  factor  equal  to  the  number  of 
processors  used. 

3.  Comparison  of  mathematical  structures  for  the  dynamic  model 

The  reduced-order  model  as  described  in  the  previous  section  is  fully  discretized:  in  time,  in  the 
state,  and  in  the  inputs.  This  mode!  can  then  be  used  to  predict  the  evolution  of  the  pair 
correlation  function  for  time-varying  but  piecewise  constant  process  inputs.  However,  one  can 
also  think  about  interpolating  between  these  discrete  points,  which  potentially  could  reduce  the 
error  due  to  the  discretization. 

A  first  method  for  interpolation  was  explored  using  the  framework  of  cel!  maps.  A  simple  cell 
map  is  a  discretized  version  of  a  continuous-state  dynamic  model,  in  which  each  discretized  state 
is  mapped  to  exactly  one  other  state  at  a  future  time  (under  a  specified  input).  This  is  the  case  in 
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our  reduced-order  model.  However  generalized  cell  maps  can  be  constructed  in  which  a  single 
state  maps  to  multiple  states  with  differing  probabilities  [4].  We  explored  this  approach  as  a 
method  to  interpolate  between  discrete  states,  but  due  to  the  coarse  discretization  of  our  state 
space,  spreading  due  to  numerical  diffusion  was  large,  causing  the  overall  prediction  error  to 
increase  in  comparison  to  simple  cell  mapping.  Generalized  cell  mapping  could  provide  an 
improvement  under  a  finer  discretization  of  the  state  space,  but  due  to  the  significant 
computation  associated  with  running  each  short  atomic-scale  simulation,  we  recommend  simple 
cell  mapping  over  generalized  cell  mapping  given  the  current  computational  technology. 

The  problem  with  the  generalized  cell  map  is  that  the  number  of  states  with  non-zero  probability 
grows  large  as  the  number  of  time  steps  advances.  To  avoid  this  problem,  while  still 
interpolating,  we  investigated  the  k-nearest  neighbors  interpolation  method  on  the  discrete  space. 
In  particular,  we  focused  on  k  =  2,  so  that  the  two  best  matching  units  of  the  current  state  of  the 
reduced-order  model  are  identified.  The  probabilities  of  each  of  these  two  states  are  computed 
based  on  the  Euclidean  distance  between  the  current  state  and  two  best  matching  states.  The  two 
best  matching  states  are  each  propagated  forward  by  one  step,  and  the  new  state  is  the  weighted 
average  of  the  two  image  cells.  In  our  investigation  of  this  approach,  we  found  that  in  most 
cases  similar  performance  was  achieved  with  simple  cell  mapping  (k  =  1)  and  k  =  2,  but  that 
occasionally  the  k  =  2  method  selected  a  second  best  matching  unit  that  was  not  a  good  match, 
and  that  caused  the  model  error  to  grow.  Thus,  for  our  relatively  coarse  discretization  of  the 
state  space,  we  saw  no  clear  advantage  to  using  k  >  I . 

4.  Quantification  of  model  reduction  error 

We  demonstrated  our  model  reduction  approach  in  a  lattice  kinetic  Monte  Carlo  simulation  of 
gallium  arsenide  deposition,  as  developed  by  Itoh  and  co-workers  [5],  Our  study  included  model 
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reduction,  open-loop  control  via  dynamic  programming,  and  a  quantification  of  the  errors  in  the 
reduced-order  models.  This  study  has  been  submitted  for  publication  in  Automatica  and  is 
included  in  this  report  as  Appendix  A. 

The  time  for  the  deposition  was  divided  into  10  intervals,  and  the  average  modeling  error  grew 
as  the  simulations  progressed,  as  would  be  expected  due  to  propagation  of  error  from  previous 
steps.  We  also  constructed  distribution  plots  for  those  errors  over  hundreds  of  trajectories,  and 
found  that  the  error  distributions  are  unimodal  and  approximately  Gaussian  in  shape  at  each  time 
step.  At  the  end  of  the  10  steps,  the  error  was  similar  in  magnitude  to  the  error  associated  with 
the  discretization  of  the  slate  space,  While  this  is  an  acceptable  and  expected  level  of  error,  we 
anticipate  that  it  will  continue  to  grow  as  the  number  of  steps  grows  larger  and  this  level  of  error 
would  become  unacceptable.  While  10  time  steps  may  seem  small  from  the  perspective  of 
simulation  and  prediction,  it  does  enable  the  consideration  of  trajectory  planning  problems  with 
10  stages  and  input  levels,  which  greatly  increases  the  number  of  possible  processing  recipes  and 
potential  material  structures.  For  example,  for  our  8  input  levels,  there  are  8UI  possible 
trajectories  that  can  be  considered.  Typically,  only  constant  inputs  are  used,  with  only  8  input 
trajectories  possible.  Furthermore,  exact  models  over  long  prediction  horizons  are  not  needed 
w'hen  feedback  control  is  used. 

Our  efforts  to  quantify  and  predict  the  error  associated  with  the  model  reduction  suggested  a  new 
approach  using  spatial  statistics  [6],  which  is  used  in  the  geology  community  under  the  name  of 
kriging.  A  new  project  has  recently  been  funded  by  AFOSR  (PI:  Gallivan)  to  develop  this 
method  in  the  context  of  dynamic  models.  In  the  past  it  has  been  used  extensively  but  for  static 
models  in  geology  and  in  mechanical  design.  It  is  expected  to  have  wide  applicability  in  error 
quantification  for  empirical  dynamics  models  including  Markov  chains  and  tabulation  models.  It 
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provides  a  unifying  framework  to  a  range  of  ad-hoc  approaches  currently  used  in  modeling  of  a 
range  of  Air  Force  applications  including  molecular  interactions,  combustion,  and  turbulent  flow. 

5.  Control 

The  majority  of  the  work  in  this  program  was  on  the  model  reduction  approach,  since  the  simple 
state-affine  structure  of  the  final  reduced-order  model  enables  the  application  of  most  existing 
methods  for  dynamic  optimization  and  control.  The  discrete  nature  of  our  state  space  enabled 
the  implementation  of  dynamic  programming  for  optimal  trajectory  generation  of  the  flux  profile 
in  the  gallium  arsenide  model.  The  time-optimal  trajectory  is  reported  in  Appendix  A,  We  also 
note  that  this  control  approach  could  be  implemented  in  real-time  using  dynamic  programming, 
since  we  can  easily  compute  off-line  the  time-optimal  trajectory  for  each  discrete  state  in  the 
state  space. 

The  implementation  of  feedback  control  requires  that  the  state  can  be  measured  or  accurately 
estimated.  Estimation  and  observers  were  not  a  primary  focus  of  the  work  supported  by  this 
grant,  although  we  did  conduct  a  study  demonstrating  that  the  reduced-order  Markov  chain 
model  can  be  used  efficiently  with  an  extended  Kalman  filter,  in  a  simulation  study  (Gallivan  in 
2005,  Computers  &  Chemical  Engineering).  Continued  w'ork  in  our  group  includes  a  joint 
modeling  and  experimental  study  for  estimation  of  microstructure  using  limited  optical  sensors 
in  a  chemical  vapor  deposition  process.  This  work  is  supported  under  an  NSF  grant  [7]. 

6.  Modeling  of  heterostructure  field  effect  transistors 

Professor  Gallivan  spent  the  summer  of  2006  working  at  Wright-Patterson  Air  Force  Base,  under 
the  support  of  a  US  Air  Force  Summer  Faculty  Fellowship.  The  goal  of  the  project  was  to 
develop  a  mathematical  model  for  the  degradation  dynamics  of  heterostructure  field  effect 


12 


transistors  (HFET)  made  from  gallium  nitride.  Based  on  the  current  understanding  of  hydrogen 
diffusion  in  gallium  nitride,  a  mechanism  for  degradation  of  AIGaN/GaN  high  electron  mobility 
transistors  was  proposed  and  simulated.  Hydrogen  is  introduced  into  the  devices  during  their 
processing  and  fabrication,  and  it  interacts  with  defects  to  electrically  passify  them.  In 
particular,  H+  diffuses  into  the  near  surface  region  of  the  AlGaN  and  GaN  during  plasma 
processing.  However,  H'  is  the  preferred  state  of  the  hydrogen  interstitial  in  undoped  GaN,  and 
thus  there  is  little  diffusion  of  hydrogen  during  device  operation.  The  traps  depassify  thermally 
during  operation,  which  is  potentially  also  enhanced  by  hot  electrons  and  electromigration.  This 
mechanism  would  not  be  reversed  by  heating,  but  if  H+  is  present,  then  heating  the  device 
uniformly  to  200-300  C  would  induce  recovery  of  a  degraded  device.  Consideration  of  the 
spatial  distributions  in  the  HFET  is  the  key  to  understanding,  modeling,  and  predicting 
degradation.  Modeling  of  the  spatial  distributions  enables  prediction  of  power  law  degradation 
dynamics  that  are  observed  in  experiments  but  have  not  previously  been  understood 
theoretically.  Experimental  validation  of  this  mechanism  was  needed,  and  experiments  at  AFRL 
are  now  in  progress.  Continued  interaction  with  AFRL  is  planned  under  the  new  grant  (PI: 
Gallivan). 

7.  Modeling  of  hyperb ranched  polymers 

Structure  development  in  highly  branched  segmented  polyurethaneureas  has  been  investigated  by 
experimental  studies  and  kinetic  Monte-Carlo  simulations.  The  experiments  were  used  to 
develop  a  minimal  set  of  reaction  kinetics  for  the  molecular  simulations.  Five  descriptors  were 
used  to  describe  the  state  of  each  polymer  (the  number  of  the  two  types  of  monomer  in  the 
molecule,  the  number  of  each  type  of  unreacted  end  groups,  and  the  number  of  intramolecular 
reactions).  Additionally  the  number  of  molecules  of  each  type  is  tracked  in  the  simulations.  The 
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simulations  have  been  used  to  suggest  new  process  trajectories,  by  changing  the  feed  profile  for 
monomer  addition,  the  dilution  of  the  system,  and  the  amount  of  monofunctional  agents  to  be 
added.  Two  papers  have  been  written,  one  of  which  is  published  (Unal  et  al.  in  2005,  Polymer) 
and  another  which  has  been  submitted  and  is  included  as  Appendix  B  of  this  report.  This  work 
has  leveraged  funding  from  the  Army  Research  Office  through  a  collaboration  with  Professor 
Tim  Long  at  Virginia  Tech.  Work  has  focused  on  validating  the  model,  so  the  model  reduction 
framework  has  not  yet  been  applied  in  this  system.  However,  part  of  this  ongoing  collaboration 
includes  the  identification  of  a  minimal  state  representation  for  molecular  architecture  in 
hyperbranched  polymers,  which  is  now  an  unsolved  problem.  The  hyperbranched  modeling 
work  has  impacted  the  model  reduction  methodology  in  a  more  indirect  way,  by  providing  an 
alternative,  and  very  different  example,  for  which  the  model  reduction  methodology  must  also  be 
applicable. 

III.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  conclusion,  the  most  significant  accomplishment  of  this  research  program  was  the 
development  of  an  automated  model  reduction  procedure  for  atomic-scale  many  body 
simulations.  The  reduction  was  demonstrated  in  a  simulation  of  gallium  arsenide  surface 
processing,  and  enabled  the  computation  of  a  time-optimal  process  trajectory  to  reach  a  desired 
surface  structure.  The  computation  required  for  generation  of  the  model  and  for  the  open-loop 
dynamic  optimization  was  tractable,  and  enabled  a  dynamic  optimization  that  would  not  have 
been  possible  using  the  full  stochastic  simulations.  Additional  materials  processing  applications 
considered  were  heterostructure  field  effect  transistors  and  hyperbranched  polymers,  and  these 
applications  also  guided  the  development  of  our  general  model  reduction  framework.  This 


14 


framework  should  have  wide  applicability  in  a  range  of  nanostructured  materials  processing 

problems. 
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Optimization  of  a  thin  film  deposition  process  using  a  dynamic 
model  extracted  from  molecular  simulations  * 
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A  bst  ract 

This  study  presents  and  demonstrates  an  algorithm  for  computing  a  dynamic  model  for  a  thin  film  deposition  process. 
The  proposed  algorithm  is  used  on  high  dimensional  Kinetic  Monte  Carlo  (KMC)  simulations  and  consists  of  applying 
principal  component  analysis  (PC A)  for  reducing  the  state  dimension,  self  organizing  map  (SOM)  for  grouping  similar  surface 
configurations  and  simple  cell  mapping  (SCM)  for  identifying  the  transitions  between  different  surface  configuration  groups. 
The  error  associated  with  this  model  reduction  approach  is  characterized  by  running  more  than  1000  test  simulations  with 
highly  dynamic  and  random  input  profiles.  The  global  error,  which  is  the  normalized  Euclidean  distance  between  the  simulated 
and  predicted  states,  is  found  out  to  be  0,006  on  average  for  the  test  simulations.  This  indicates  that  our  reduced  order 
dynamic  model,  which  was  developed  using  a  rather  small  simulation  set,  was  able  to  accurately  predict  the  evolution  of  the 
film  mic restructure  for  much  larger  simulation  sets  and  a  wide  range  of  process  conditions.  Minimization  of  the  deposition 
time  to  reach  a  desired  film  structure  has  also  been  achieved  using  this  model.  Hence,  our  study  showed  that  the  proposed 
algorithm  is  useful  for  extracting  dynamic  models  from  high  dimensional  and  noisy  molecular  simulation  data. 


Key  words:  Dynamic  Modelling;  Model  reduction;  Modelling  errors;  Optimal  control;  Order  reduction. 


1  Introduction 

Thin  film  deposition  is  a  critical  step  in  manufacturing 
integrated  circuits  and  MEMS  devices.  As  device  size 
gets  smaller,  film  thickness  and  tolerance  approach  the 
atomic  scale,  so  high  quality  thin  films  with  uniform  and 
smooth  surfaces  arc  desired  for  high  device  performance. 
Surface  processing  is  commonly  used  in  developing  inte¬ 
grated  circuits  in  order  to  build  features  with  dimensions 
of  100  ran  and  below  [1,2].  Some  other  applications  of 
surface  processing  arc  mechanical  coatings  [3],  thermal 
coatings  [4]  and  MEMS  devices  [5].  The  structure  of  a 
thin  film  is  usually  a  a  strong  function  of  process  inputs 
such  as  pressure,  plasma  power  and  temperature.  There¬ 
fore,  process  models  that  describe  the  relationship  be¬ 
tween  inputs  and  outputs  are  needed  to  control  the  tni- 
crostructure  of  thin  films.  Since  continuum  assumptions 
arc  not  valid  at  such  small  scales,  developing  accurate 
continuum  models  from  material  and  energy  balances 
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is  not  possible.  An  alternative  approach  is  simulating 
the  dynamics  of  the  system  by  means  of  molecular  sim¬ 
ulations,  such  as  molecular  dynamics  (MD)  and  Monte 
Carlo  (MC)  simulations.  In  MD  simulations,  Newton's 
equations  of  motion  are  solved  for  the  position  of  each 
atom  in  a  system.  Given  the  initial  location  of  each  atom, 
a  potential  energy  function  is  used  to  compute  the  inter¬ 
action  between  atom  pairs  and  this  information  allows 
the  computation  of  the  trajectory  of  each  atom  over  a 
time  interval.  On  the  other  hand,  Monte  Carlo  (MC) 
simulation  is  a  stochastic  method  which  generates  dif¬ 
ferent  discrete  configurations  of  a  system  by  randomly 
changing  the  position,  orientation  and  conformation  of 
the  atoms  or  molecules  in  the  system. 

The  specific  application  considered  in  this  paper  is  the 
epitaxial  growth  of  gallium  arsenide  (GaAs),  GaAs  is 
generally  grown  by  ultra-high  vacuum  molecular  beam 
epitaxy  (MBE).  There  arc  many  advantages  of  using 
GaAs  over  silicon.  It  has  a  higher  saturated  electron  ve¬ 
locity  and  higher  electron  mobility,  allowing  devices  to 
function  at  frequencies  excess  of  250  GHz.  Also,  GaAs 
devices  generate  less  noise  than  silicon  devices.  Having 
a  direct  band  gap,  GaAs  can  be  used  to  emit  light  unlike 
silicon,  which  has  an  indirect  band  gap  and  is  a  poor  light 
emitter.  MBE  deposition  of  GaAs  occurs  with  a  rate  of 
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approximately  one  layer  per  second.  Hence,  the  growth 
morphology  in  this  process  develops  in  the  order  of  sec¬ 
onds.  Even  with  a  further  increase  in  computer  power, 
such  a  time  scale  will  not  be  accessible  to  simulations 
using  MD.  Even  though  MD  can  capture  the  atomic  vi¬ 
brations  in  the  order  of  picoseconds,  slower  events  like 
an  atom  overcoming  an  energy  barrier  and  moving  to  a 
new  site  in  the  crystal  lattice  will  be  infrequent  when 
this  technique  is  used.  However,  the  slower  events  will  be 
dominant  in  terms  of  capturing  the  dynamics  of  the  pro¬ 
cess,  Therefore,  we  use  Kinetic  Monte  Carlo  (KMC)  sim¬ 
ulations  to  reach  the  time  scales  of  these  slower  events. 

MC  simulations  are  typically  used  to  compute  equi¬ 
librium  properties*  In  contrast,  Kinetic  Monte  Carlo 
(KMC)  simulations  are  also  able  to  describe  evolution 
in  time.  The  KMC  algorithm  was  first  proposed  by 
Bortz  et  al.  [6]  and  then  transformed  by  Maksyrn  [7]  to 
be  used  in  MBE  growth.  It  is  based  on  the  assumption 
that  the  surface  evolution  is  the  result  of  a  scries  of  dis¬ 
crete  events.  These  events  arc  classified  as  adsorption 
(atom  attaching  to  the  surface),  surface  diffusion  (atom 
hopping)  and  desorption  (atom  leaving  the  surface).  In 
a  KMC  simulation,  the  initial  arrangement  of  atoms  is 
specified  and  the  possible  transitions  from  this  configu¬ 
ration  are  evaluated.  Then,  based  on  generated  random 
numbers,  one  transition  and  Its  location  is  selected. 
After  the  transition  is  executed,  the  simulation  clock 
is  incremented  using  another  random  number.  A  KMC 
simulation  model  on  a  zincblcndc  crystal  structure  has 
been  developed  to  describe  epitaxial  growth  of  GaAs 
[8]*  This  model  includes  the  rates  (derived  from  exper¬ 
iments)  of  over  one  thousand  possible  events  taking 
place  on  the  surface. 

Several  approaches  have  been  developed  in  the  past  few 
years  to  integrate  molecular  simulations  with  dynamic 
analysis  and  optimization*  One  approach  is  the  deriva¬ 
tion  of  stochastic  time  evolution  equations  from  the 
probabilistic  master  equation  (11],  under  assumptions 
that  are  mostly  applicable  for  well-mixed  reacting  sys¬ 
tems.  An  alternative  approach  for  constructing  stochas¬ 
tic  differential  equations  from  molecular  simulations  of 
film  growth  is  the  generation  of  stochastic  spatially  dis¬ 
tributed  PDEs  for  the  height  profile  of  a  surface  [12]* 
These  models  describe  continuum  behavior  and  dynam¬ 
ics,  and  thus  are  not  appropriate  for  modeling  surface 
structure  at  atomic  scales*  This  modeling  approach, 
with  stochastic  PDEs,  has  been  used  recently  to  control 
the  surface  roughness  in  a  thin  film  deposition  process 
on  a  one-dimensional  lattice  [13]*  On  the  other  hand, 
controller  design  based  on  molecular  simulations  is  an 
active  area  of  research  [14].  In  a  recent  study  [15],  a  re¬ 
duced  order  stochastic  model  obtained  from  KMC  and 
finite  difference  simulation  data  has  been  used  to  design 
a  feedback- feed  forward  controller  in  order  to  maintain 
tlic  current  density  during  the  copper  electrodeposition 
process  at  a  constant  level. 


In  order  to  explicitly  model  and  predict  discrete  atomic 
scale  structure,  a  model  is  required  that  docs  not  av¬ 
erage  over  the  small  length  scales.  One  approach  to 
address  this  issue  is  equation-free  computing*  which 
was  first  used  for  stability  and  bifurcation  analysis  [9]* 
A  low-dimensional  system  state  was  assumed  (based 
on  macroscopic  arguments),  and  short  simulations  were 
run  from  specific  initial  conditions  to  approximate  time- 
derivatives*  More  recently,  this  method  has  also  been 
applied  for  optimization  [10],  However,  the  reduction 
in  computational  time  achieved  by  this  method,  com¬ 
pared  to  running  full  molecular  simulations,  may  not  be 
sufficient  to  make  the  approach  practical  when  many 
predictions  over  long  times  intervals  arc  required* 

In  an  alternative  approach  presented  by  Gal U van  and 
Murray  [16],  molecular  simulations  were  used  to  con¬ 
struct  Markov  models,  with  discrete  states  describing 
groups  of  similar  configurations.  This  grouping  strategy 
was  based  on  the  similarity  of  the  roughness  values*  and 
enabled  computation  of  the  optimal  temperature  pro¬ 
file  by  penalizing  the  surface  roughness  and  temperature 
changes  during  the  thin  film  deposition.  The  construc¬ 
tion  of  this  explicit  low- order  model  reduced  the  com¬ 
putational  load  by  four  orders  of  magnitude,  when  com¬ 
pared  with  full  molecular  simulations.  This  significant 
reduction  was  critical  for  making  the  dynamic  simula¬ 
tion  and  optimization  feasible  over  macroscopic  process¬ 
ing  times. 

The  construction  of  a  dynamic  model  relies  on  the  exis¬ 
tence  and  knowledge  of  tbc  system  state*  While  a  low- 
order  state  was  selected  in  previous  studies  using  physi¬ 
cal  arguments  [9,16],  the  detailed  structure  in  a  molecu¬ 
lar  simulation  may  not  even  have  a  low-order  representa¬ 
tion,  Additionally,  small  structural  features  may  have  a 
large  effect  on  the  time  evolution  of  the  system.  In  order 
to  address  this  issue,  a  previous  study  by  Oguz  and  Gal- 
livan  [17]  proposed  the  use  of  high-dimensional  step-step 
correlation  functions,  which  provided  a  more  detailed 
state  description  of  the  film  surface,  compared  to  the 
surface  roughness  alone.  This  study  demonstrated  the 
implementation  of  principal  component  analysis  (PCA) 
for  reducing  the  state  dimension  and  self  organizing  map 
(SOM)  for  automated  grouping  of  the  similar  states  in 
the  state  space.  In  a  more  recent  study,  Varshncy  mid 
Armaou  [18]  also  used  spatial  correlation  functions  to 
characterize  the  state  of  their  film  growth  simulation, 
and  used  equation- free  computing  to  simulate  the  dy¬ 
namics  [9]* 

The  next  section  explains  our  modeling  approach,  which 
consists  of  state  reduction  and  model  identification. 
State  reduction  is  achieved  by  reducing  the  dimension 
of  the  simulation  data  using  principal  component  anal¬ 
ysis  (PCA)  and  grouping  similar  surface  configurations 
using  a  seif  organizing  map  (SOM).  After  establishing 
a  discrete  state  space,  transitions  between  different  sur¬ 
face  configuration  groups  have  been  identified  by  simple 
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Fig.  I .  Schematic  of  the  modeling  approach 

cell  mapping  (SCM)  for  model  identification  purposes. 
In  Section  3,  we  evaluate  the  performance  of  the  process 
model  by  computing  the  local  and  global  prediction 
error  under  highly  dynamic  input  profiles.  Section  4 
describes  the  minimization  of  the  film  deposition  time 
to  reach  an  optimal  film  structure.  Finally,  Section  5 
provides  the  conclusions  of  this  study. 

2  Modeling  Approach 

Our  modeling  approach  t  which  consists  of  four  steps T 
is  illustrated  in  Figure  1 .  Characterization  of  the  state 
space,  reduction  of  state  dimension,  discretization  of  the 
state  space,  and  model  identification  gives  us  a  reduced 
order  dynamic  model. 

2.  t  Characterization  of  the  state  space 

A  step-step  correlation  (SSC)  function  is  used  to  de¬ 
scribe  the  microscopic  state  of  the  system  during  the 
simulations.  This  function  gives  the  distance  and  orien¬ 
tation  between  pairs  of  steps  on  the  surface,  where  a  step 
is  defined  as  a  change  in  height  from  one  atomic  surface 
site  to  the  next.  Since  the  occurrence  of  every  type  of 
step  pair  (up- up,  up- down,  down-up  and  down- down)  is 
counted  on  the  surface  in  both  directions  on  the  Carte¬ 
sian  coordinate  axes,  the  spatial  correlation  function  is 
high-dimensional  and  may  contain  redundant  informa¬ 
tion.  The  data  is  also  noisy.  Noise  is  reduced  by  perform¬ 
ing  multiple  realizations  under  identical  conditions  and 
averaging  the  results.  After  averaging,  PC  A  is  used  to 
determine  the  number  of  independent  variables  needed 
to  fully  determine  the  SSC  function.  This  technique  is 
widely  used  to  eliminate  linear  correlations  among  the 
variables  in  data  sets.  It  is  a  crucial  step  in  our  study. 
The  reduced  state  dimension  makes  it  possible  to  con¬ 
struct  a  compact  dynamic  model. 

The  KMC  simulations  have  been  carried  out  using  the 
kinetic  barriers  calculated  by  It  oh  jS]  with  the  following 
parameters: 

*  Growth  temperature:  580  C 

*  Film  deposition  interval:  0.20  monolayers  (ML) 


■  Incident  arsenic  dimer  (AS2)  fiux:  0.4  ML/s 

*  Incident  gallium  (Ga)  flux:  Varied  between  0.06-0.20 
ML/s  (0.06,  0.08,  0.10,...,  0.20)  flux  range  where  the 
model  is  valid. 

•  Lattice  size:  300x300  (90000  surface  atoms)  with  pe¬ 
riodic  boundary  conditions 

The  starting  surface  in  our  simulations  is  the  thermo¬ 
dynamic  ground  state  of  GaAs(001),  the  /?2(2  x  4)  re¬ 
construction  which  prevails  in  a  wide  range  of  growth 
conditions. 

Snapshots  of  surfaces  have  been  recorded  at  surface  in¬ 
crements  of  0.01  ML,  starting  from  an  initial  surface  in 
the  ;32(2  x  4)  configuration,  up  to  0.20  ML  coverage. 
Eight  constant  input  simulations  were  performed  at  Ga 
fluxes  of  0.06  ML/s,  0.08  ML/s,  0.10  ML/s,...,  0.20  ML/s 
and  we  will  refer  this  as  Training  Simulation  Set,  1.  In 
Training  Simulation  Set  2,  we  again  have  8  simulations, 
but  this  time  the  flux  is  kept  constant  up  to  0.10  ML  cov¬ 
erage  (middle  of  the  deposition),  and  the  flux  is  shifted 
to  a  different  value  at  that  coverage  point.  To  explore  the 
KMC  state  space  even  further,  we  perforin  60  additional 
simulations  (Third  Training  Simulation  Set),  where  two 
flux  shifts  arc  made  at  0.07  and  0. 14  ML  coverage  points. 

The  SSC  function  has  been  computet]  as  S  €  R'*5  for 
each  snapshot.  S  is  obtained  by  combining  the  16  dif¬ 
ferent  portions  of  the  SSC  function,  which  has  been 
evaluated  for  four  different  types  of  step  pairs  in  four 
different  directions  on  the  surface  (4x4—16  combina¬ 
tions).  Each  portion  has  300  variables,  which  is  equal 
to  the  lattice  size.  Therefore,  ds  has  a  value  of  4800 
(300  x  16  —  4800).  We  have  a  total  number  of  snap¬ 
shots  na=1521  (76  simulations,  20  snapshots  for  each 
simulation,  plus  the  snapshot  of  the  initial  state  of  the 
system).  S  of  all  snapshots  in  the  training  data  arc  col¬ 
lected  in  D  £  R4800*1531.  Before  performing  PCA,  D  is 
first  transformed  into  D *  €  r180D*1521  with  the  follow¬ 
ing  elements: 


~Dfncan,t.  (1) 

Here,  Dme( m  €  R4800  is  defined  as: 


(SjLi  Du) 

Dmean.i  =  - L  (i  =  1, 2,  ..,1800).  (2) 

Tlj, 


In  order  to  complete  the  pre-processing  of  D,  Df  is  trans¬ 
formed  into  D "  €  R4800*1  . 
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Dfj  =  Dynam  (i  =  1,2,  ,.4800,  j  =  1,2,  ..1521) 

(3) 


D«td  6  K^00  is  defined  as: 


Hence,  by  pre-processing  D,  the  variance  of  each  variable 
in  the  SSC  function  is  made  equally  important  before 
PC  A  is  carried  out.  In  other  words,  small  features  with 
low  variance  will  not  bo  neglected. 


2.2  Reduction  of  the  state  dimension 

Principal  component  analysis  (PC A)  has  been  per¬ 
formed  through  the  singular  %'aluc  decomposition  of 
matrix  Dn  by  computing  the  singular  values  of  this 
matrix.  The  squares  of  the  singular  values  correspond 
to  the  eigenvalues  of  the  covariance  matrix  of  Df\  and 
the  ratio  of  each  eigenvalue  to  the  sum  of  all  eigen¬ 
values  (normalized  eigenvalue)  is  plotted  against  the 
principal  components.  The  point  on  the  plot,  where  a 
sudden  decay  of  the  normalized  eigenvalue  is  seen,  gives 
the  minimum  number  of  principal  components  (the 
minimum  dimension  n)  that  can  reconstruct  the  data 
effectively.  At  this  point,  we  project  each  snapshot’s 
SSC  function  onto  the  first  n  principal  components  Df\ 
We  define  x  €  Rn  its  the  coefficient  set  with  these  new 
coordinates.  Each  x  characterizes  a  particular  snapshot. 
As  a  result  of  PCA,  small  features  in  the  SSC  function, 
which  do  not  contribute  to  the  variance  of  the  Dn  signif¬ 
icantly,  are  eliminated.  Wc  note  that  this  could  create  a 
problem  later  while  grouping  similar  surface  structures 
according  to  first  n  principal  components,  because  these 
small  features  may  possess  valuable  information  about 
the  differences  between  the  surface  structures  and  their 
evolution  in  time. 


2.3  Discretization  of  state  space 

After  performing  PCA  and  obtaining  a  coefficient  set  to 
characterize  each  surface  snapshot,  we  use  self  organizing 
map  (SOM)  to  eliminate  some  of  the  nonlinear  correla¬ 
tions  within  the  data  set  by  employing  MATLAB’s  SOM 
Toolbox.  Snapshots  with  very  similar  microstructurcs 
are  grouped  by  SOM,  and  snapshots  in  the  same  group 
are  viewed  as  equivalent  when  identifying  the  dynamic 
model  in  the  next  step.  The  PCA  data  is  discretized  as  a 
part  of  this  step,  since  snapshots  are  grouped  among  the 
nodes  of  SOM.  As  a  result,  the  number  of  discrete  states 
is  finite.  The  computational  load  for  the  model  identi¬ 
fication  is  also  reduced  with  the  grouping  achieved  by 


SOM.  Because,  during  model  identification,  rather  than 
the  transitions  between  each  surface  structure,  only  the 
transitions  between  the  groups  are  computed. 

Various  training  procedures  for  SOM  arc  described  else¬ 
where  [19].  SOM  training  is  accomplished  by  sequen¬ 
tial  or  batch-wise  comparison  of  snapshots  to  the  overall 
map,  and  then  updating  the  prototype  vectors  v  e  RTi 
of  the  map  nodes  to  match  and  organize  the  snapshots. 
Once  the  map  is  trained,  each  snapshot  is  associated 
with  the  node  that  has  the  closest  prototype  vector,  as 
measured  by  the  Euclidean  distance.  The  default  num¬ 
ber  of  nodes  in  the  map  is  computed  by  the  SOM  tool¬ 
box  using  a  heuristic  formula  which  is  a  function  of  the 
number  of  rows  (snapshots)  in  the  data  matrix.  The  map 
size  can  also  be  increased  to  provide  a  finer  discretiza¬ 
tion  of  the  state  space.  Another  important  SOM  param¬ 
eter,  the  ratio  of  the  sidclcngths  of  the  map,  is  set  equal 
to  the  ratio  of  the  two  largest  principal  components  of 
the  data  matrix, 

2.3  A  Results 

In  this  study,  PCA  is  used  to  find  the  minimum  di¬ 
mension  that  can  represent  the  microscopic  state  of  the 
surface  snapshots  recorded  during  the  simulations.  This 
method  consists  of  computing  the  variance  captured  by 
each  principal  component  of  the  entire  data  set  and  se¬ 
lecting  the  most  important  principal  components  that 
can  reconstruct  the  data  well.  Figure  2  shows  the  nor¬ 
malized  eigenvalues  of  each  principal  component.  A  knee 
shape  is  observed  after  the  second  principal  component 
and  the  first  twTo  principal  components  capture  nearly 
all  of  the  variance  within  the  data  set. 


However,  data  reconstructions  with  2  and  even  more 
principal  components  showed  us  that  at  least  5  com¬ 
ponents  arc  needed  in  order  to  effectively  represent  our 
data.  Specifically,  5  principal  components  reconstructed 
the  small  clusters  of  atoms  on  the  surface  (with  a  size  of 
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less  than  20  lattice  units)  much  better  than  2  principal 
components.  Figure  3  illustrates  this  comparison. 


Fig.  3.  Comparison  of  the  reconstructions  using  2  and  5 
modes  with  the  original  data  for  snapshot  61.  This  surface 
is  obtained  under  0, 10  Mb/s  Ga  flux  and  has  a  0,20  ML 
surface  coverage.  The  region  enclosed  with  the  gray  dashed 
line  represents  the  surface  structures  with  a  size  of  less  than 
20  lattice  units. 

As  already  defined  in  Section  2.2,  the  coefficient  set  x 
for  each  surface  snapshot  is  obtained  by  taking  the  inner 
product  of  the  snapshot’s  SSC  function  with  the  n  =  5 
principal  components  of  the  data  matrix  Dtf .  This  re¬ 
duced  represent  at  ion*  which  includes  the  coefficient  sets 
of  every  snapshot,  is  used  to  train  an  SOM.  For  train¬ 
ing.  we  used  the  default  map  size,  which  is  based  on  the 
number  of  surface  snapshots  in  the  training  data.  Each 
SOM  node  lias  its  own  prototype  vector  having  the  same 
dimension  as  each  snapshot  in  the  data  set  (a  =  5),  As 
a  result,  each  map  node  also  has  its  own  SSC  function 
S,  which  is  the  PC  A  reconstruction  of  its  prototype  vec¬ 
tor  v.  Thus,  each  map  node  is  associated  with  a  type  of 
surface. 

The  quality  of  the  resulting  map  can  be  evaluated  by 
calculating  the  average  quantization  error  over  the  input 
data,  given  by  ||x  -  v0My||2,  where  vbmu  is  the  proto¬ 
type  vector  for  the  best  matching  unit,  for  the  snapshot 
with  the  projection  vector  x  [19].  In  words,  the  Eucli- 
adcan  distance  between  each  surface  snapshot  and  its 
best  matching  SOM  node  is  computed,  and  the  average 
of  these  distances  gives  the  average  quantization  error. 

In  order  to  see  the  effects  of  training  data  on  the  map 
quality,  we  generated  multiple  maps  while  keeping  the 
map  size  at  192  nodes.  For  each  map,  the  training  phase 
has  taken  approximately  3  minutes  using  a  Pentium  4 
processor  with  a  speed  of  3  GHz.  The  first  map  SOMl  is 
trained  with  161  surface  snapshots  (from  Training  Sim¬ 
ulation  Set  1),  the  second  map  SO  M2  is  trained  with 
321  surface  snapshots  (from  Training  Simulation  Sets  1 
and  2)  and  SOM3  is  trained  using  all  of  the  1521  sur¬ 
face  snapshots  (from  Training  Simulation  Sets  1 ,  2  and 


3),  Tabic  1  shows  the  statistics  of  three  maps  with  192 
nodes,  which  were  trained  using  these  three  different 
data  sets.  SOMl  and  SOM2  have  similar  statistics,  ex¬ 
cept  that  SOM2,  which  has  been  trained  with  a  larger 
data  set,  has  a  larger  number  of  surface  snapshot  groups. 
This  suggests  that  Training  Simulation  Set  2t  where  the 
flux  was  changed  in  the  middle  of  the  simulations  at  0. 10 
ML  surface  coverage,  enabled  us  to  sec  surface  struc¬ 
tures  that  had  not  been  accessible  through  constant  in¬ 
put  simulations  (Training  Simulation  Set  1).  On  Table 
1,  as  we  move  from  SOM2  to  SOM3,  wc  again  observe 
an  increase  in  the  number  of  surface  snapshot  groups. 
Therefore,  we  conclude  that  Training  Simulation  Set  3 
(during  which  the  flux  has  been  changed  twice  during  the 
simulations)  explored  the  KMC  state  space  even  further, 
producing  additional  surface  structures  compared  with 
Training  Simulation  Sets  1  and  2.  SO  M3  has  an  average 
quantization  error  of  5, 1 T  whereas  the  average  prototype 
vector  size  is  53.0.  If  needed,  the  average  quantization 
error  can  be  reduced  by  adding  more  nodes  to  the  map, 
but  we  were  able  to  obtain  a  useful  model  with  SO  M3, 

Another  SOM  metric  is  the  topographic  error,  which  is 
the  proportion  of  all  the  data  vectors  for  which  first  and 
second  best  matching  units  arc  not  adjacent  on  the  map 
[19].  SOM3  also  has  a  topographic  error  of  0,01,  which 
means  that  for  only  15  of  the  1521  snapshots,  the  first 
and  second  best  matching  units  were  not  adjacent  on  the 
SOM.  Topology  preservation  is  not  required  for  a  good 
state  representation ,  but  it  aids  in  visualization  of  the 
dynamics  on  the  map. 

At  this  point,  three  SOMs  have  been  generated.  Only 
one  of  them,  SOM3t  will  be  used  to  identify  a  dynamic 
model,  since  it  has  been  trained  with  the  most,  extensive 
amount  of  training  data. 

Tabic  1 

Statistics  of  three  SOMs  trained  using  different  data  sets 


Map 

Training 

sim,  sets 

Quant. 

error 

Froto. 

vcc.  size 

Topog. 

error 

Number 

of  groups 

SOMl 

1 

4,1 

54.1 

0.01 

119 

SO  M2 

1,2 

4,8 

54.4 

0,02 

142 

SO  M3 

1,2,3 

5,1 

53.0 

0,01 

175 

2-4  Model  Identification 
2-4-  /  Simple  Cell  Mapping 

Characterization  of  the  surface  snapshots  from  all  of 
the  training  simulations  gives  the  system’s  state  space. 
As  explained  in  the  previous  section,  the  surface  snap¬ 
shots  arc  grouped  into  discrete  states  using  SOM,  ac¬ 
cording  to  their  similarities  in  terms  of  surface  morphol¬ 
ogy.  These  groups  represent  the  cells,  and  simple  cell 
mapping  (SCM)  or  generalized  cell  mapping  (GCM)  can 
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be  used  for  obtaining  a  global  view  of  the  behavior  of 
the  system  [21],  A  cell  map  is  formed  by  dividing  the 
state  space  into  a  finite  number  of  cells  (using  SOM) 
and  approximating  the  behavior  of  the  system  by  means 
of  transitions  between  the  cells.  SCM  is  a  determinis¬ 
tic  approach  in  which  each  cell  is  mapped  into  exactly 
one  other  cell  for  a  particular  input.  In  GCM,  each  cell 
can  have  several  image  cells.  In  other  words,  a  cell  can 
be  mapped  to  several  other  cells.  Using  this  stochastic 
approach,  we  can  assign  a  probability  of  the  system  be¬ 
ing  in  a  cell  at  a  specific  time  and  extract  the  dynamic 
properties.  In  the  current  study,  wc  implement  the  SCM 
approach  to  identify  the  dynamic  model  for  the  pro  cess , 
because  SCM  provides  a  deterministic  model  to  describe 
the  evolution  of  the  surface  structure  under  different  ma¬ 
terial  dux  profiles.  The  following  arc  the  steps  used  for 
model  identification: 

*  One  surface  configuration  is  selected  from  a  particular 
map  node. 

*  A  new  simulation  is  started  from  that  surface  config¬ 
uration  and  run  for  an  incremental  coverage  interval 
with  one  of  the  flux  settings. 

*  The  final  surface  structure  is  recorded  and  the  config¬ 
uration  group,  which  has  the  closest  microstructure 
to  that  final  structure,  is  obtained.  In  other  words, 
the  flux-dependent  transition  (from  one  configuration 
group  to  another)  is  found. 

When  the  procedure  described  above  is  repeated  for  ev¬ 
ery  map  node  and  each  flux  setting,  wc  identify  a  reduced 
order  dynamic  model.  The  identification  procedure  re¬ 
quires  running  146  x  8  =  1168  KMC  simulations  (for 
146  map  nodes  and  8  flux  settings),  which  takes  approx¬ 
imately  7.3  days  using  a  computer  duster  made  up  of  16 
computers,  each  having  an  Intel  Xcon  processor  with  a 
speed  of  2.66  GHz. 

The  reduced  order  model  is  described  as: 


p\j  +  1]  =  -A(u[j])p[j]  (5) 

where  j  is  the  time  step  (or  coverage  level}  number, 
u|j]  6  R  is  the  input  (or  flux)  value,  p  €  Rm  is  the  proba¬ 
bility  vector  describing  the  configuration  group  that  the 
system  is  currently  in  and  A  €  Rmxm  is  the  transition 
matrix.  The  value  of  m  (dimension  of  p)  is  equal  to  the 
number  of  configuration  groups  in  the  state  space  (the 
number  of  cells  in  the  cell  map).  We  define  p  as: 

I  if  i  =  /; 

0  if  i  ^  L 

Here,  l  is  the  number  of  the  current  system  configura¬ 
tion.  The  transition  matrix  A  is  a  function  of  the  input. 


Therefore,  in  the  reduced  order  model,  we  have  eight 
transition  matrices  For  eight  flux  settings.  The  sum  of 
the  dements  in  each  column  of  A  is  equal  to  one,  bo 
cause  for  each  column,  only  one  clement  of  A  is  non-zero, 
whereas  the  other  dements  are  equal  to  zero.  In  other 
words,  there  is  only  one  possible  transition  from  each 
configuration  group  under  a  particular  input.  This  con¬ 
struction  of  A  is  due  to  the  deterministic  nature  of  the 
simple  cell  map.  If  the  system  makes  a  transition  from 
configuration  s  to  t  at  time  step  j  with  the  flux  setting 
u[j],  this  transition  can  be  represented  by  a  transition 
matrix  with  the  following  properties: 

1  if  k  =  t; 

0  if  k  #  f. 

At  time  step  j ,  the  surface  properties  of  the  system  are 
given  by  x\j]  (set  of  projection  coefficients  defined  in 
Section  2.2): 


4j]=*p[j]  (6) 

where  X  €  Rn*m  has  the  v  €  Rn  (prototype 
vector)  of  each  configuration  group  (SOM  node). 
In  other  words,  v  of  the  configuration  group  i  is 
v*  —  { A  * ht i  j  Xj hi,  Xst, }. 

2,4*2  K -nearest  neighbor  algorithm 

Cell  mapping  can  be  generalized  to  enable  Interpolation 
between  cells.  In  order  to  carry  out  the  interpolation, 
wc  use  the  k-ncarcst  neighbor  algorithm  [20]  and  predict 
the  evolution  of  x  for  a  given  flux  profile.  When  k=L 
this  algorithm  produces  results  identical  to  simple  cell 
mapping  and  no  interpolation  is  made.  The  following 
steps  represent  the  k- nearest  neighbor  algorithm  when 
k— 2: 

(1)  Let  xf,fEj  be  the  vector  representing  the  surface  state 
at  an  initial  film  coverage.  This  vector  is  compared 
with  all  the  prototype  vectors  on  the  SOM  and  its 
best  matching  unit  (BMU)  and  second  BMU  arc 
found  based  on  the  Euclidean  distances  between 
these  vectors.  Let  these  these  two  SOM  nodes  have 
prototype  vectors  v  i  and  v<2 ,  where  v  x  and  Vg  repre¬ 
sent  the  first  BMU  and  second  BMU,  respectively. 

(2)  Compute  the  distances  d\  =  V||[a  and  ds  = 

IxoM-Vflla- 

(3)  Compute  the  weights  associated  with  these  two  dis¬ 
tance  values  ziq  =  4-  l/d?)  and  — 

(iM)/(i/dt  +  i  M)- 

(4)  Predict  the  system  state  at  the  next  step  from 
xnew  —  uq  .v3  +  u^  v-t,  where  is  the  state  vec¬ 
tor  prediction  for  the  next  step,  v3  and  are  the 
prototype  vectors  of  the  nodes  determined  from 
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SCM  (vt  transitions  into  v3,  and  v2  transitions 
into  v4  under  the  particular  flux  setting). 

(5)  Set  x0u  =  xnew  and  repeat  steps  2-4  to  make  pre¬ 
dictions  at  higher  film  coverage  values.  At  each  cov¬ 
erage  value,  the  predicted  x  can  bo  converted  to  the 
high  dimensional  SSC  function  by  multiplying  the 
elements  of  x  with  the  eigenvectors  obtained  from 
the  training  data  D 

The  algorithm  described  above  is  the  k-ncarcst  neigh¬ 
bor  algorithm  with  k=2.  Here,  the  interpolation  between 
prototype  vectors  provide  the  opportunity  to  improve 
the  accuracy  of  the  prediction.  On  the  other  hand,  when 
k=l.  k- nearest  neighbor  algorithm  is  identical  to  SCM. 
In  that,  case,  the  second  BMU  is  not  identified  and  the 
weights  are  not  computed. 

3  Performance  Evaluation  of  the  Dynamic 
Model 

3 * l  Local  Error  Quantification 

The  core  hypothesis  of  our  model  reduction  approach 
is  that  the  surface  structures  in  the  same  configuration 
group  (SOM  node)  should  show  similar  dynamic  behav¬ 
ior  under  identical  process  input  (material  flux).  SCM, 
which  is  used  to  extract  a  dynamic  process  model  in  this 
study,  is  also  useful  for  testing  this  hypothesis.  We  define 
the  cell  mapping  error  (CM E)  to  quantify  the  different 
dynamic  behaviors  of  surface  structures  which  belong  to 
the  same  configuration  group.  The  following  is  the  pro¬ 
cedure  we  used  to  And  CM  E\ 

(1)  Randomly,  select  one  surface  structure  from  a  map 
node  and  identify  where  this  structure  is  mapped 
on  the  SOM  under  a  particular  flux  setting  (first 
cell  mapping). 

(2)  Randomly,  select  another  surface  structure  from 
the  same  map  node  and  run  another  simulation 
starting  with  this  new  surface  structure  and  per¬ 
form  SCM  again  (second  cell  mapping) . 

(3)  Compute  CME  =  |jSj  -  S2|h/(|[Si  +  S3||3/2)f 
where  Si  and  S2  arc  the  SSC  functions  of  the  SOM 
nodes  coming  from  the  first  and  second  cell  map¬ 
ping,  respectively*  It  should  be  noted  that  these 
functions  are  reconstructed  from  the  prototype  vec¬ 
tors  of  the  map  nodes, 

(4)  Repeat  steps  1-3  for  all  SOM  nodes  under  all  flux 
settings. 

CM E  is  computed  by  picking  two  different  surface  struc¬ 
tures  from  each  node  and  performing  SCM  under  each 
of  the  eight  flux  settings.  In  order  to  characterize  the 
distribution  of  CME,  its  cumulative  distribution  func¬ 
tion  is  computed.  This  function  is  defined  as  CDF(y)  = 
P(CME  <  y)  and  gives  the  value  of  the  probability 
that  CME  <  y  for  a  given  y.  Figure  4  shows  the  CDF 
of  the  CAfE .  52%  of  the  mappings  are  identical,  with 


both  surfaces  evolving  to  the  same  map  node.  In  those 
cases,  CME  -  0.  Also,  with  a0.f)0  probability  CME  < 
0.0075,  supporting  the  fact  that  there  is  a  very  strong 
chance  for  surface  structures  in  the  same  group  to  show 
similar  dynamic  behavior.  In  order  to  reduce  CM  E *  a 
larger  SOM  can  be  used.  The  tradeoff  is  that  this  would 
increase  the  dimension  of  the  cel)  map  and  the  compu¬ 
tational  load  associated  with  the  system  identification 
step. 


Fig,  4.  Cumulative  distribution  function  of  the  cell  mapping 
error 

3,2  Global  Error  Quantification 

A  dynamic  model  is  constructed  by  finding  the  flux  de¬ 
pendent  transitions  between  the  nodes  of  SOM3  under 
all  flux  settings.  In  order  to  test  the  predictive  ability  of 
our  model,  1210  test  simulations  are  performed.  In  these 
simulations,  we  split  the  film  coverage  domain  (from  0  to 
0.20  ML)  into  10  equal  intervals.  The  gallium  flux  value 
is  kept  constant  in  each  interval.  However,  its  value  ran¬ 
domly  changes  when  moving  from  one  coverage  interval 
to  the  next.  This  strategy  provides  very'  different  flux 
profiles  than  the  ones  used  to  generate  the  training  data, 
which  involved  a  maximum  of  two  flux  switches.  Figure 
5  shows  the  flux  profile  for  one  of  these  simulations  dur¬ 
ing  which  the  flux  is  randomly  changed  at  each  step. 

As  described  in  Section  2T  surface  snapshots  from  the 
test  simulations  at  different  coverage  levels  (10  surface 
snapshots  from  0.02  ML  to  0.20  ML  film  coverage)  arc 
characterized  using  SSC  functions,  and  the  coefficient 
set  x  €  of  cadi  surface  snapshot  is  computed.  Then* 
each  x  is  matched  with  an  SOM  node  based  on  the  cri¬ 
teria  of  minimum  Euclidean  distance  between  x  and  the 
prototype  vectors  of  the  SOM  nodes.  In  other  words, 
the  best  matching  unit  for  each  x  is  sought.  For  the  test 
simulation  with  the  flux  profile  shown  in  Figure  5t  the 
predicted  trajectory  and  the  KMC  simulation  trajectory 
arc  given  in  Figure  6.  Here,  the  prediction  comes  from 
the  k- nearest  neighbor  algorithm  (with  k=I)  and  each 
hexagon  represents  an  SOM  node  corresponding  to  a 
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Fig.  6.  Trajectories  of  the  KMC  test  simulation  (red  line) 
and  the  prediction  (black  line)  on  SO  M3 

surface  structure  group.  The  simulation  follows  a  path 
starting  from  the  SOM  Node  1  (initial  surface  structure) 
and  moves  to  the  right  hand  side  of  the  SOM  as  the 
film  coverage  increases.  Since  each  simulation  starts  from 
SOM  Node  1  and  the  transitions  from  each  node  under 
all  flux  settings  arc  known,  the  prediction  of  the  simu¬ 
lation  trajectory  (using  the  dynamic  model)  is  straight¬ 
forward.  Figure  6  indicates  that  there  is  an  agreement 
between  the  predicted  trajectory  and  the  trajectory  of 
the  KMC  simulation.  Similar  results  were  obtained  for 
all  test  simulations.  In  some  instances,  these  trajectories 
passed  through  neighboring  SOM  nodes.  However,  as 
Table  1  indicates.  SO  M3  has  a  very  low  topographic  er¬ 
ror.  Hence,  neighboring  SOM  nodes  are  similar,  so  slight 
differences  between  these  trajectories  do  not  jeopardize 
the  accuracy  of  the  state  predictions. 

SOM  Node  186  represents  the  predicted  film  structure 
at  the  final  coverage  value  (0.20  ML).  Figure  7  shows  the 
SSC  function  reconstructed  from  the  prototype  vector 
of  this  SOM  node.  The  reconstruction  agrees  with  the 
KMC  simulation’s  SSC  function.  Hence,  our  dynamic 
model  is  capable  of  predicting  the  final  film  structure 
quite  well  This  figure  also  illustrates  that  the  noise  in 
the  KMC  simulation  data  is  considerably  reduced  when 
the  SSC  function  is  reconstructed  from  the  prototype 
vector. 

In  order  to  evaluate  the  performance  of  the  reduced  order 
dynamic  model,  wc  quantify  the  prediction  error  for  the 


Fig.  7.  Reconstruction  of  the  SSC  function  with  the  proto¬ 
type  vector  of  the  SOM  Node  186  and  the  original  KMC 
simulation  data  at  final  film  coverage 

tost  simulations,  which  were  not  included  in  the  training 
data.  The  global  (multi-step)  error  associated  with  the 
predicted  SSC  function  is  defined  as: 


Esse  =  l|S,  -  Splla/IIS.Ia,  (7) 

where  Sa  is  computed  from  the  KMC  simulation  data 
and  S;)  is  reconstructed  from  the  predicted  state  vector 
v  (prototype  vector  of  the  SOM  node).  The  values  of 
ESse  at  different  film  coverage  levels  (for  the  test  sim¬ 
ulation  with  the  flux  profile  in  Figure  5)  arc  less  than 
0.018  according  to  Figure  8.  This  figure  compares  the 
error  in  the  predictions  made  using  the  k-ncarcst  neigh¬ 
bor  algorithm  with  different  k  values.  For  this  particular 
test  simulation,  predictions  arc  less  accurate  with  k— 2. 
However,  for  more  than  90%  of  the  test  simulations,  pre¬ 
dictions  did  not  change  when  k  value  was  increased  from 
1  to  2.  This  is  because,  at  low  film  coverage  levels,  differ¬ 
ent  surface  configuration  groups  (corresponding  to  dif¬ 
ferent  SOM  nodes)  do  not  possess  significantly  different 
surface  features  (or  different  x)  and  they  are  mapped 
to  the  same  SOM  node  at  subsequent  coverage  levels. 
Therefore,  throughout  the  rest  of  this  section,  we  only 
report  the  results  obtained  with  k— 1,  which  is  identical 
to  simple  cell  mapping. 

As  a  part  of  the  global  error  quantification,  we  also  com¬ 
puted  Esse  at  nine  evenly  distributed  film  coverage  val¬ 
ues  (0.04  ML,  0.06  ML,...,  0.20  ML)  for  three  large  test 
simulation  sets  (with  random  input  profiles),  none  of 
which  had  been  in  the  SOM  training  data.  It  should  be 
noted  that  the  SOM  node  1,  which  represents  the  initial 
film  structure  at  0  ML.  maps  to  the  same  SOM  node 
at  0.02  ML  with  ail  flux  settings.  Hence,  the  prediction 
error  at  0.02  ML  is  not  computed.  The  three  test  sim¬ 
ulation  sets  include  300,  600  and  1210  simulations,  re¬ 
spectively.  Figure  9  compares  the  mean  values  of  Esse 
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Fig.  8.  The  evolution  of  Esse  for  the  test  simulation  with 
k=l  and  k— 2 


Fig.  9.  The  mean  value  of  Esse  at  different  film  coverage 
levels  for  three  test  simulation  sets  and  the  training  data 

at  different  film  coverage  levels  for  these  three  test  sim¬ 
ulation  sets  and  also  the  Training  Simulation  Set  3.  The 
mean  Esse  is  below  0.012  at  all  coverage  levels.  It  is 
interesting  to  note  that  the  dynamic  model  was  able  to 
predict  the  final  film  structure  in  the  training  simula¬ 
tions  with  less  error  than  the  test  simulations.  This  is 
because  of  the  fact  that  the  data  from  the  training  sim¬ 
ulations  were  used  as  an  input  in  the  system  identifica¬ 
tion,  where  as  our  dynamic  model  was  not  as  familiar 
with  some  of  the  surface  structures  produced  during  the 
test  simulations,  especially  at  high  film  coverage* 

As  shown  in  Table  2,  the  mean  values  of  Esse  for  the 
predictions  associated  with  the  test  simulation  sets  are 
around  0*006,  but  the  standard  deviation  values  arc  com¬ 
parable  to  the  mean  values  indicating  a  wide  distribu¬ 
tion  of  the  prediction  error.  In  order  to  get  a  more  clear 
idea  about  the  distribution  of  E$sc,  its  cumulative  dis¬ 
tribution  function  (CDF)  is  computed  For  the  three  test 
simulation  sets  and  also  the  Training  Simulation  Set  3* 
Figure  10  shows  that  the  CDF  curves  of  the  test  sim- 


Fig.  10.  Cumulative  distribution  function  of  the  Esse  for 
three  test  simulation  sets 

illation  sets  are  very  similar  and  enlarging  the  size  of 
the  simulation  set  does  not  change  the  distribution  of 
the  prediction  error  significantly.  Also,  the  probability 
of  having  an  Esse  less  than  0.01  is  around  90%,  which 
again  supports  the  fact  that  the  reduced  order  dynamic 
model  has  a  good  prediction  capability. 

Table  2 

Mean  and  standard  deviation  values  of  Esse  for  three  test 
simulation  sets 


Test  simulation  set 

Mean 

Standard  deviation 

1 

0.0057 

0.0035 

2 

0,0058 

0*0037 

3 

0.0058 

0,0037 

In  the  last  part  of  this  section,  we  compare  different  types 
of  prediction  error  to  understand  the  major  factors  con¬ 
tributing  to  the  inaccuracies  in  the  model  predictions. 
We  define  two  alternatives  to  Esse  in  order  to  investi¬ 
gate  the  effect  of  normalization  (using  Dmcail  and  DMtd 
defined  in  Section  2.1.)  on  the  prediction  error.  The  first 
alternative  to  Esse  is 


Esse*  =  IIS',  -S'p||a/||SMa  {8} 

where  S'*  and  S'P  arc  the  normalized  versions  of  S9  and 

V 

S'a,i  =  S»,i  -  Dmtan,i  ( i  -  1, 2,  ...4800)  (9) 


S'p,i  =  SPil  -  Dmea„'i  (i  =  1,2,  ...4800)  (10) 

is  already  defined  in  Section  2*1*  The  second 
alternative  to  Esse  is: 
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EW  =  l|S",-S%||a/||S"J2  (11) 

whore  S"fl  and  S"p  are  the  normalized  versions  of  S'* 
and  S'p,  respectively: 


S".,i  =  s 

(i  =  1,2,  ...4800) 

(12) 

S"p,i  =  S'p,i/D.td,i 

(i  =  1,2,  ...4800) 

(13) 

jDatcti  was  also  defined  in  Section  2.1.  Table  3  shows  that 
mean  value  of  Esse "  (0.2675)  is  much  higher  than  that 
of  Egsc *  (0.1758).  Comparison  of  the  CDF  of  Esse '  and 
Esse”  given  in  Figure  11  also  shows  a  higher  Esse 
for  any  value  of  the  CDF.  This  indicates  that  features 
with  large  variance  are  predicted  more  accurately  than 
the  ones  with  small  variance,  It  follows  then,  that  the 
features  with  small  variance  arc  more  noise- corrupted  in 
the  KMC  simulations. 

The  quantization  error,  which  is  defined  in  Section  2.3.1, 
provides  a  minimum  bound  on  the  prediction  error.  The 
average  quantization  error  (£^)  for  SGM3  is  0. 1094,  This 
error  is  due  to  the  discretization  of  the  state  space,  and 
docs  not  include  any  additional  error  propagated  from 
one  step  to  the  next  one  through  the  dynamics  of  the 
thin  film  deposition  process.  Therefore,  we  define  an¬ 
other  kind  of  prediction  error  as 


E*  =  l|x*  -xp||2/]|xj2.  04) 

1  Icrc,  xp  €  R5  represents  the  predicted  state  and  xa  € 
R5  is  computed  from  the  KMC  simulation  data: 


xPi,  =  S"p  ■  U,r 

(<-1,2,. ..5) 

(15) 

x,,  =  S",  -  U,T 

(i  =  1,2,  ...5) 

(16) 

where  U,  e  R4800  is  the  i*h  principal  component  of  Dtf. 

As  shown  in  Table  3,  the  mean  value  of  Ex,  which  is  com¬ 
puted  from  the  prediction  of  states  in  Test  Simulation 
Set  3,  is  twice  as  high  as  E v  incurred  during  the  training 
of  SO  M3.  From  this  result,  it  can  be  concluded  that  the 
propagation  error  due  to  dynamics  is  not  negligible  and 
accounts  for  approximately  half  of  Ex *  Hence,  the  other 
half  of  the  prediction  error  conies  from  the  discretization 
of  the  state  space.  Figure  1 1 ,  which  compares  the  CDF 
of  Ex  and  E<v  also  supports  this  conclusion.  Here,  for 
any  particular  value  of  error,  the  probability  of  having 
that  particular  value  of  Ex  is  higher  than  that  of  E,r 


Tabic  3 

Mean  values  of  Eq,  ET,  Es$cf  and  Esse"-  which  is  com¬ 
puted  for  each  data  vector  and  the  prototype  vector  of  the 
data  vector's  best  matching  unit  on  SOM,  is  the  difference 
between  data  vector  and  prototype  vector  during  SOM  train¬ 
ing*  The  other  errors  arc  the  normalized  versions  of  Esse- 


Error  type 

Mean 

0.1 0(M 

Ex 

0,2277 

Esse 

0.0057 

Esse* 

0.1 758 

ESSc" 

0.2675 

Fig.  1 1 .  Cumulative  distribution  function  of  different  types 
of  error 

4  Optimization  of  the  final  film  structure  and 
the  deposition  time 

In  this  section,  our  dynamic  model  is  used  to  minimize 
the  deposition  time  to  reach  a  particular  surface  config¬ 
uration.  One  desired  structure  is  a  very  regular  surface 
with  many  large  islands  and  a  very  low'  Ga  adatom  (iso¬ 
lated  Ga  atom)  density.  This  surface  structure  can  he 
identified  by  minimizing: 


F  =  Oi-  bid  (17) 

where  a,i,  bi,  and  c;  arc  the  values  of  Ga  adatom  den- 
sitv,  typical  island  size  and  the  number  of  islands  with 
the  typical  island  size  for  surface  configuration  i,  respec¬ 
tively,  From  the  Training  Simulation  Sets  1 ,  2  and  3,  we 
extract  a*,  6*,  and  ct  values  of  the  surfaces  at  0.20  ML 
surface  coverage. 

According  to  equation  (17),  the  optimal  surface  struc¬ 
ture  is  Snapshot  861  (from  Simulation  43),  Snapshot  861 
is  matched  with  the  Node  182  of  SOM3  during  the  train¬ 
ing.  This  node  is  accessible  through  a  constant  input 
simulation  with  a  flux  of  0.08  ML/s.  Identification  of  the 
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flux  profile  that  would  reach  Node  182  in  the  minimum 
amount  of  time  can  be  posed  as  adynamic  programming 
problem  [22]. 

Let  the  decision  variables  dq(q  =  1,  2, 3,  ...11)  be  the  im¬ 
mediate  destinations  on  different  stages.  In  this  prob¬ 
lem,  we  have  1 1  stages  for  11  film  coverage  levels  (0  ML, 
0.02  ML,  0.04  ML,,., 0,20  ML)  and  dq  corresponds  to  the 
map  node  number  at  stage  q.  Thus  the  route  (trajec¬ 
tory)  of  the  deposition  is  di,  dn,  where  dj  =  1 

and  du  —  182  since  the  initial  surface  structure  is  rep¬ 
resented  by  SOM  Node  1  and  the  final  (optimal)  surface 
structure  is  in  Node  182. 

Let  fq(s.  dq)  be  the  total  cost  of  the  best  overall  policy 
for  the  remaining  stages,  given  that  we  arc  in  state  s 
(number  of  the  map  node  we  arc  currently  in),  ready  to 
start  stage  q ,  and  dq  (the  number  of  the  map  node  we 
are  moving  to)  is  selected  as  our  immediate  destination. 
Hero,  the  total  cost  is  the  deposition  time  and  each  de¬ 
position  interval  is  0,02  ML  long.  Given  s  and  q T  let  d* 
denote  any  value  of  dq  that  minimizes  fq(stdq)  and  let 
/*(s)  be  the  corresponding  minimum  value  of  /^(s,  d,;). 
Thus, 


/'(s)  =  min  fq(s,d9)  =  fq(s,  d’)  (18) 

where 

Msydq)  =c,'d,  +  /,+i(s,d,+1)  (19) 

Here,  the  cost  ca^i  is  the  time  incurred  while  moving 
from  S  to  dq,  given  by: 


c9,49  —  LfF9tdqi  (20) 

where  L  is  the  length  of  a  single  coverage  interval  (0.02 
ML)  and  F3  <irf  is  the  value  of  the  gallium  flux  that  pro¬ 
vides  a  transition  from  state  s  to  dcr  The  ultimate  desti¬ 
nation  reached  at  the  end  of  stage  11,  fy{  (182)  =  0,  The 
objective  is  to  find  fi(l)  and  the  corresponding  route. 
Dynamic  programming  can  solve  this  problem  by  succcs- 
sivdy  finding /9*(s),  fs(s)-fi(s)  and  using  /£( s) 
to  solve  for  fi(s).  This  is  achieved  by  eliminating  some 
of  the  suboptimal  trajectories  as  wo  move  from  /*n(s) 
to  Because  of  the  limited  state  space  obtained  by 

grouping  similar  surface  con  figuration  groups,  it  is  possi¬ 
ble  to  solve  this  dynamic  optimization  problem  (finding 
the  optimal  flux  profile)  using  exhaustive  enumeration 
(without  eliminating  the  suboptimal  paths  at  each  step) 
in  a  short  amount  of  time. 

Wc  used  eight  flux  settings  (0.06 ,0.08 ,...0.20  ML/s)  and 
found  the  optimal  input  profile  that  would  give  us  min¬ 
imum  deposition  time  to  get  to  SOM  Node  182.  Each 
input  profile  is  a  sequence  of  10  flux  values,  and  there 
is  a  flux  value  for  each  coverage  interval.  Having  eight 
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Fig.  12.  Optimal  flux  profile  computed  by  the  dynamic  model 

flux  settings  and  10  coverage  intervals,  there  are  810  pos¬ 
sible  input  profiles.  According  to  the  dynamic  model, 
only  20%  of  these  profiles  arc  able  to  reach  SOM  Node 
182.  Running  each  simulation  takes  about  24  hrs  utiliz¬ 
ing  an  Intel  Xcon  processor  with  a  speed  of  2. 66  GHz, 
so  it  would  have  taken  2.9  million  years  to  run  all  of  the 
810  simulations  with  a  single  processor.  However,  using 
our  dynamic  model,  it  took  only  5  minutes  to  predict 
the  evolution  of  the  film  structure  during  these  simula¬ 
tions.  The  minimum  cost  was  obtained  with  the  input 
profile  shown  in  Figure  12,  Figure  13  shows  the  real  and 
the  estimated  trajectories  of  this  KMC  simulation  on 
SO  M3.  Again,  wfc  have  a  good  agreement  between  these 
two  trajectories.  This  particular  input  profile  provided 
a  48%  reduction  in  the  deposition  time  to  reach  opti¬ 
mal  structure  when  compared  with  the  constant  input 
KMC  simulation  under  0.08  ML/s  Ga  flux.  The  values 
of  Esse*  Es sc*  and  Esse "  for  the  prediction  of  the  fi¬ 
nal  film  structure  are  0.0159,  0.2207  and  0.231 2t  respec¬ 
tively.  These  values  are  within  the  range  of  cumulative 
distribution  functions  shown  in  Figures  10  and  11.  In  or¬ 
der  to  visualize  the  accuracy  of  this  prediction,  we  also 
plotted  a  portion  of  the  SSC  functions,  which  belonged 
to  the  actual  KMC  simulation  data,  its  best  matching 
unit  on  the  map  (Node  184)  and  Node  182  (predicted 
film  structure).  According  to  Figure  14,  the  simulation 
data  is  much  noisier  than  the  reconstructions.  Also,  the 
lower  peak  value  of  the  SSC  function  (around  radius=  12) 
of  the  SOM  Node  182  is  slightly  off  compared  to  the  one 
coming  from  the  KMC  simulation.  However,  the  plateau 
corresponding  to  the  number  of  step  pairs  which  are  dis¬ 
tanced  with  more  than  20  lattice  units  is  captured  well 
with  the  prediction.  Hence,  these  two  SSC  functions  are 
very  similar.  This  similarity  indicates  that  the  dynamic 
model,  once  again,  docs  a  good  job  in  terms  of  predict* 
ing  the  final  film  structure.  Furthermore,  this  prediction 
can  be  improved  by  increasing  the  size  of  the  training 
data  set  and  number  of  SOM  nodes  in  the  cell  map,  or 
possibly  by  placing  a  greater  weight  on  the  important, 
but  small  sized  features,  during  the  training  of  the  SOM, 

5  Conclusions 

This  study  describes  an  approach  to  understanding  and 
modeling  process  dynamics  from  molecular  simulations. 
The  proposed  algorithm  is  used  on  high  dimensional  Ki¬ 
netic  Monte  Carlo  (KMC)  simulations  of  epitaxial  GaAs 
thin  film  deposition.  This  algorithm  consists  of  three 
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Fig.  13,  Trajectories  of  the  KMC  simulation  (red  line)  with 
the  optimal  flux  profile  and  the  prediction  (black  line)  on 
SO  M3 


Fig.  14.  Reconstructions  of  the  SSC  function  with  the  pro¬ 
totype  vectors  of  the  SOM  nodes  182,  184  and  the  original 
simulation  data 

steps:  applying  principal  component  analysis  (PCA),  self 
organizing  map  (SOM)  and  simple  cell  mapping  (SCM) 
to  identify  a  dynamic  process  model.  First,  a  spatial 
correlation  function  is  used  to  describe  the  state  space 
of  the  system  by  the  characterization  of  surface  snap¬ 
shots  generated  under  different  material  flux  profiles. 
Then,  a  minimum  state  dimension  that  can  represent  the 
system  state  is  found  using  PC  A.  After  this  reduction 
stepT  SOM  is  employed  to  group  similar  surface  struc¬ 
tures.  SOM  results  led  to  the  ident  ification  of  a  dynamic 
model  through  the  computation  of  flux- dependent  tran¬ 
sitions  between  surface  configuration  groups  with  simple 
cell  mapping  (SCM),  Analysis  of  the  cell  mapping  error 
(CME),  which  is  the  one-step  prediction  error,  showed 
that  the  structures  within  the  same  configuration  groups 
show  very  similar  dynamic  behavior  under  same  input 
conditions.  The  global  error  associated  with  this  model 
reduction  approach  has  also  been  characterized  using 


1210  test  simulations  with  highly  dynamic  input  profiles 
and  turned  out  to  be  fairly  low  (less  than  0.006  on  aver¬ 
age  for  10890  predictions).  Furthermore,  the  minimiza¬ 
tion  of  the  deposition  time  to  reach  a  desired  film  struc¬ 
ture  has  also  been  achieved  using  the  compact  dynamic 
model.  Our  study  shows  that  the  proposed  algorithm  is 
useful  for  extracting  reduced  order  dynamic  models  from 
high  dimensional  and  noisy  molecular  simulation  data. 
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Interpretation  of  molecular  structure  and  kinetics  in  melt  condensation  of  A; 
oligomers,  Bj  monomers,  and  monofunctional  reagents 

C.  Oguz,  S.  Unal,  T.  E.  Long,  and  M.  A.  Gallivan 

Abstract 

This  manuscript  applies  kinetic  Monte  Carlo  simulations  to  interpret  experimental 
measurements,  in  the  polymerization  of  hyperbranched  poly(ether  esters)s  in  a  melt 
condensation  of  A  2  oligomers  and  B3  monomers.  Building  on  the  analytical  modeling  of 
Flory  and  Stockmayer,  additional  effects  of  cycle  formation,  unequal  reactivities,  and 
end-capping  reagents  are  added  into  the  simulations,  to  describe  A2+B3  polymerization  in 
the  absence  of  a  solvent.  The  experimental  data  has  been  published  separately1,  and  here 
it  is  compared  to  the  model  predictions  in  order  to  quantitatively  assess  whether  the  data 
is  consistent  with  these  models.  Based  on  the  modeling,  we  conclude  that  cycle 
formation  is  negligible,  suppression  of  the  third  B-group  is  insignificant,  and  the  mobility 
of  the  free  B3  monomer  leads  to  enhancement  of  its  reaction  rate.  The  addition  of  the 
monofunctional  end-capping  reagents  does  not  necessarily  lead  to  suppression  of 
branching  in  the  A2+B3  system,  and  depends  sensitively  on  the  stoichiometry  of  the 
reactants. 


1/26 


Appendix  B 

Submitted  to  Macromolecules 

1.  Introduction 

The  shape  and  topology  of  organic  molecules  has  a  profound  effect  on  their  properties.2 
In  the  last  two  decades,  synthetic  polymer  chemists  have  introduced  a  new  class  of  highly 
branched  macromolecules  which  are  composed  of  multifunctional  monomers,  and  are 
classified  as  either  dendrimers  or  hyperbranched  polymers,  Dendrimers  are  typically 
synthesized  using  multi-step  reactions  and  they  offer  superior  control  of  molecular  size, 
shape,  and  functionality.  On  the  other  hand,  hyperbranched  polymers  are  less  ordered  but 
are  easier  to  synthesize. 

Discussion  of  the  synthetic  methodologies  for  preparation  of  a  wide  range  of 
hyperbranched  and  dendritic  polymers  can  be  found  in  several  extensive  review  articles.3' 
5  Many  of  the  past  experimental  studies  have  focused  on  the  synthesis  and 
characterization  of  ABn  type  monomers  (n>l),  particularly  with  AB2  monomers.3"6 
However,  very  few  of  the  AB„  type  monomers  are  commercially  available  due  to  their 
lack  of  symmetrical  functionality  and  tendency  to  react  prematurely.4  As  a  result.  A2+B3 
polymerization  has  recently  been  the  subject  of  extensive  research7'14  since  it  provides  an 
alternative  and  more  convenient  way  to  synthesize  highly  branched  polymers.  In  contrast 
to  polymerization  of  ABn  type  monomers,  these  systems  offer  a  wider  range  of  molecular 
structures  depending  on  the  monomeric  types  and  processing  conditions.  For  example, 
A2+B3  polymerization  has  been  performed  by  heating  a  mixture  of  Aj+Bj,1  and  also  by 
drop-wise  addition  of  A2  into  Bj  IS.  The  molar  ratio  of  A 2  to  B3  can  also  be  varied  4-  l6,  l7. 
Moreover,  our  recent  efforts  have  demonstrated  the  opportunity  to  control  the  distance 
between  branch  points  through  the  judicious  selection  of  various  telechelic  oligomers  as 
A2. 
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Modeling  studies  have  long  been  used  to  explain  experimental  observations  in  an 
adhoc  fashion,  and  modeling  is  often  used  to  steer  the  discovery  of  synthetic  methods  and 
the  formation  of  novel  architectures.  Early  work  of  Flory16  and  Stockmayer1'3,  |l}  on  the 
step  growth  of  multifunctional  monomers  was  based  on  the  assumption  that  there  was  no 
cycle  formation  in  these  polymerization  processes,  which  enabled  the  calculation  of 
molecular  weight  using  an  infinite  series  solution.  The  models  of  Flory  and  Stockmayer 
are  useful  because  they  enable  quantitative  predictions  of  polymer  properties.  While  it  is 
intuitive  that  a  branched  polymer  that  is  composed  of  A2  and  Bj  monomers  will 
ultimately  gel,  the  model  generates  the  quantitative  prediction  of  gelation  at  87%  A 
conversion  for  the  stoichiometry  of  Aj:!^  =  1:1. 19  More  recently,  this  method  has  been 
extended  to  include  the  effect  of  cycle  forming  reactions  on  the  gel  point.211 

In  order  to  predict  the  time-evolution  of  polymerization,  kinetic  models  based  on 
mass-action  kinetics  are  often  employed.  Ordinary  differential  equations  are  used  to 
describe  the  concentration  of  each  type  of  branching  unit.21,  22  For  example,  a 
trifunctional  monomer  may  exist  in  four  states:  no  reactions  (free),  one  reaction 
(terminal),  two  reactions  (linear),  and  fully  reacted  (dendritic).  These  models  are 
typically  of  low  dimension  and  may  be  solved  analytically  in  some  cases,  or  numerically 
in  others.  A  disadvantage  of  this  modeling  approach  is  that  no  information  is  provided  on 
the  molecular  weight  distribution  of  the  resulting  polymer. 

When  prediction  of  the  molecular  weight  and  its  distribution  is  also  desired, 
population  balance  models  are  commonly  used,  in  which  the  concentration  of  polymers 
of  each  possible  size  is  computed  17,  2t>  23‘25.  These  models  therefore  have  high 
dimension.  For  systems  of  linear  polymers  composed  of  a  single  monomer  type,  the 
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kinetics  are  fully  described  by  the  number  of  monomers  in  the  polymer.  Thus,  if  one 
models  the  concentration  of  polymers  up  to  a  maximum  size  of  1000  monomeric  units, 
then  the  dimension  of  the  model  is  also  1000.  However,  in  branched  polymers,  more 
information  is  needed  to  describe  the  kinetics,  such  as  the  number  of  reactive  end  groups 
in  the  polymer.  When  several  descriptors  are  needed  to  describe  each  polymer,  the 
dimension  of  the  population  balance  mode!  grows  rapidly.' '  For  the  model  considered  in 
this  manuscript,  six  descriptors  are  used  to  describe  each  polymer:  the  number  of  A2 
monomers,  the  number  of  B3  monomers,  the  number  of  unreacted  A  groups  in  the 
polymer,  the  number  and  type  of  unreacted  B  groups  (linear  or  terminal),  and  the  number 
of  end-capping  agents.  Solving  a  population  balance  model  with  six  descriptors  would  be 
very  intensive  computationally. 

Generating  functions  and  the  method  of  moments  are  often  used  to  reduce  the 
dimension  of  these  population  balance  models.  However,  this  approach  becomes  much 
more  difficult  as  the  number  of  descriptors  grows,  since  moments  must  be  included  for 
each  descriptor.  The  difficulty  of  this  approach  is  illustrated  by  the  recent  paper  of 
Dusek,  Duskova-Srmckova.  and  Voit17,  in  which  unequal  reactivities  and  monofunctional 
reagents  were  considered  separately,  but  not  simultaneously. 

Cycle  formation  was  neglected  by  Flory  and  Stockmayer,  but  it  is  a  major  factor 
in  the  structural  development  of  dendritic  and  hyperbranched  polymers  23, 25':s.  The 
method  of  moments  has  been  extended  to  describe  cycle  formation  in  hyperbranched 
polymers  composed  of  AB2  monomers  and  the  AB2+B3  system  '9.  This  extension  is 
enabled  because  in  these  systems,  a  maximum  of  one  cycle  is  possible,  at  which  point  the 
polymer  cannot  grow  further.  In  contrast,  in  the  A2+B3  system,  there  is  no  limit  on  the 
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maximum  number  of  cycles.  Infinite  series  solutions  have  been  developed  to  predict  the 
gel  point  in  the  A2+B3  system,  but  do  not  include  additional  effects  such  as  unequal 
reactivities  of  the  groups  in  B3,  and  the  effect  of  monofunctional  reagents.211 

As  an  alternative  to  mass  action  kinetic  models  that  use  concentration  variables, 
Monte  Carlo  simulations  have  been  performed  for  hyperbranched  polymers,  so  that  more 
realistic  kinetics  can  be  included  and  so  that  that  structural  information  can  be  obtained. 
For  example,  the  Wiener  index  can  be  computed  for  each  polymer  in  the  simulation, 
which  is  related  to  viscosity.30  In  the  Monte  Carlo  simulations,  individual  monomers  are 
reacted  with  each  other  to  build  up  the  polymer,  using  random  numbers  to  select  each 
event.  Each  monomer  is  tracked  throughout  the  simulation,  and  information  regarding  its 
connection  to  other  monomers  is  stored.  Monte  Carlo  simulations  may  be  used  to 
describe  only  the  connectivity  of  the  polymers,  without  describing  their  spatial  positions, 
is.25,29.31  Qr  iattjce  Monte  Carlo  simulations  can  be  performed  in  which  each  monomer  is 
associated  with  a  spatial  position  in  the  lattice.32  Direct  comparison  between 
experimental  data  and  Monte  Carlo  simulations  has  been  limited  to  date,  but  these 
simulations  can  be  valuable  in  interpreting  experimental  results.15' 25  As  computational 
resources  grow,  Monte  Carlo  simulations  become  an  increasingly  attractive  alternative  to 
population  balance  modeling,  since  their  major  drawback  has  been  the  amount  of 
computation  required. 

Another  modeling  approach  is  to  describe  the  structural  development  of 
hyperbranched  polymers  using  atomistically  detailed  algorithms.1  30  These  studies  are 
limited  to  the  growth  of  single  polymers  at  a  time  because  of  the  high  computational  cost, 
so  they  are  not  as  useful  for  predicting  the  molecular  weight  distribution.  However,  if  the 
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conformation  of  the  polymer  changes  throughout  the  reaction  and  this  strongly  influences 
the  kinetics,  then  it  may  be  necessary  to  include  this  level  of  detail. 

Recent  studies  employ  syntheses  of  hyperbranched  polymers  via  A2+B3  both  in  a 
solution7'9,  12  and  also  in  the  melt  phase1  33  with  no  solvent.  Herein,  we  combine 
experimental  and  computational  efforts  to  understand  the  polymer  structural  development 
in  a  melt,  especially  branching  and  the  onset  of  gelation.  The  experimental  data  has  been 
published  previously,1  and  this  paper  employs  a  Monte  Carlo  simulation  to  interpret  the 
experimental  findings  and  to  eludicate  the  underlying  kinetics.  In  a  previous  study,  we 
used  a  simpler  kinetic  model  to  explore  the  role  of  cyclization  in  the  drop-wise  addition 
of  A2  into  Bj  in  a  solvent.15  Herein  we  consider  a  similar  statistical  framework  to  explore 
a  wider  range  of  phenomena,  for  a  batch  A2+B3  reaction  in  a  melt.  The  Monte  Carlo 
simulations  are  used  to  interpret  the  experimentally  measured  number-averaged 
molecular  weight,  weight-averaged  molecular  weight,  and  the  density  of  branched  units, 
by  considering  the  effects  of  the  cyclization  reactions,  unequal  reactivities,  and  end¬ 
capping  on  the  structure  development  of  hyperbranched  polymers. 

II.  Experimental 

Highly  branched  polyfether  ester)s  were  synthesized  in  the  melt  phase  using  an 
oligomeric  A2  +  B3  polymerization  strategy.  Condensation  of  polypropylene  glycol)  (A2 
oligomer)  and  trimethyl  1,3,5-benzenetricarboxylate  (B3  monomer)  generated  highly 
branched  structures.  The  conversion  and  degree  of  branching  were  measured  with  1 M 
NMR  spectroscopy,  and  the  molecular  weight  (number-average,  weight-average,  and 
polydispersity)  was  characterized  using  size  exclusion  chromatography  (SEC),  at  six 
points  throughout  the  polymerization  process.  Additional  experiments  were  also 
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performed  in  which  two  different  mono  functional  endcapping  reagents  were  added  for 
the  purpose  of  delaying  the  gel  point.  This  experimental  work  has  been  published 
previously,1  and  should  be  consulted  for  further  detail  on  the  experimental  procedures. 
The  present  work  interprets  these  experiments  through  a  comparison  with  Monte  Carlo 
simulations,  utilizing  alternative  assumptions  and  kinetic  models.  Because  the  error  in 
the  SEC  measurements  is  approximately  10%,  our  goal  is  not  to  achieve  exact  agreement 
between  the  experiments  and  modeling,  but  to  compare  the  trends  and  magnitudes. 

III.  Model  and  Simulations 

This  paper  focuses  on  the  use  of  a  model  to  interpret  the  experimental  data.  Each 
simulation  starts  with  N  monomers  of  A2  and  N  monomers  of  Bj  in  the  system,  since  the 
monomers  were  polymerized  with  a  1 : 1  molar  ratio  during  the  experiments.  At  each  step 
in  the  simulation,  all  of  the  available  (unreacted)  functional  groups  are  listed.  Then, 
using  a  random  number,  an  A  group  and  a  B  group  are  selected  from  the  list  and  the 
reaction  is  executed.  This  is  followed  by  updating  the  list  of  available  A  and  B  groups, 
molecular  weights  of  the  molecules,  and  the  number  of  dendritic,  linear  and  terminal 
units  in  the  system. 

The  probability  of  selecting  a  particular  pair  of  A  and  B  groups  is  proportional  to 
the  reaction  rate  for  that  pair,  which  yields  the  correct  time-evolution  of  the  system/'1 
Therefore,  a  model  is  needed  for  the  reaction  rates  of  various  events.  The  first  effect  that 
is  considered  in  this  effort  is  the  formation  of  cycles  through  intramolecular  reactions. 
Our  polymerization  was  performed  in  the  melt,  which  was  hypothesized  to  minimize 
cycle  formation  due  to  high  concentration  of  reactants,  so  cyclization  reactions  are  taken 
into  account  in  the  simulations  to  address  this  hypothesis.  While  selecting  an  A  group 
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(from  A2  oligomers)  and  a  B  group  (from  B3  monomers)  for  reaction  at  each  simulation 
step,  some  pairs  are  favored  more  than  the  others.  If  the  reaction  of  an  A-B  pair  leads  to 
cycle  formation,  the  selection  probability  of  that  particular  pair  is  promoted  by  the 
cyclization  parameter  y,  such  that  y=(kc/knc)/N,  where  kc  is  the  rate  of  cyclization 
reactions  for  each  A-B  pair  and  kne  is  the  rate  of  non-cyclization  reactions  for  each  A-B 
pair.  N  appears  in  the  parameter  y  since  the  number  of  intermolecular  reactions  in  the 
simulations  grows  as  N2,  while  intramolecular  reactions  are  initially  proportional  to  N. 
Although  y  is  constant  throughout  each  simulation,  the  overall  number  of  cycle-forming 
reactions  increases  with  conversion,  because  the  number  of  possible  cycle-forming 
reactions  increases  with  molecular  weight.  We  use  different  values  of  y  in  the 
simulations  to  explore  the  effect  of  cycle  formation  on  the  resulting  polymer. 

In  addition  to  cyclization  and  non-cyclization  reactions,  we  also  consider  end- 
capping  reactions  between  E  groups  (from  monofunctional  end-capping  reagents)  and  B 
groups,  with  a  rate  constant  of  kc.  The  ratio  e  =  ke/knc  is  then  the  second  parameter  in  our 
kinetic  model.  In  a  second  set  of  simulations,  y  and  e  have  been  varied  to  observe  their 
effects  on  the  development  of  molecular  characteristics  such  as  number  average 
molecular  weight  (M„),  weight  average  molecular  weight  (Mw),  polvdispersity  index 
(PDI  =  and  the  fraction  of  dendritic  units  (fo).  fo  is  calculated  using:  fo  -  D  / 

(D+T+L),  where  D,  L  and  T  indicate  the  number  of  dendritic,  linear  and  terminal  units  in 
the  system.  We  plot  fo  in  the  simulation  results  (as  opposed  to  fL  or  ft),  since  fo  is  the 
quantity  extracted  from  the  'H  NMR  measurements.1 

Depending  on  the  properties  of  the  monomers,  such  as  their  molar  mass  or 
electrostatic  interactions,  unequal  reactivities  of  their  groups  can  also  affect  the  structural 
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development  of  hyperbranched  polymers.  End  groups  usually  have  higher  reactivities 
than  the  groups  along  the  length  of  the  chains  because  of  the  lower  kinetic  excluded 
volume  effect.35  This  would  cause  the  linear  units  to  have  lower  reactivities  than  the 
terminal  units.  In  order  to  simulate  unequal  reactivities  of  the  B  groups,  we  define  a  third 
parameter  p.  For  each  unreacted  B  group  in  a  B3  monomer,  w-e  check  the  reaction  state 
of  the  other  two  B  groups  in  the  monomer.  We  then  consider  three  possible  cases  of 
unequal  reactivity.  For  this  purpose,  we  define  kj  to  be  the  rate  of  reaction  of  a  B-group 
in  a  free  B3  monomer,  k2  to  be  the  rate  of  a  B-group  in  a  terminal  unit,  and  k3  for  a  B- 
group  in  a  linear  unit.  We  assign  reaction  rates  in  the  ratio  of  p  =  k|/k2  =  k:/kj.  k|  is 
expected  to  be  enhanced  relative  to  kj  due  to  the  greater  mobility  of  the  free  B3  monomer 
and  its  ability  to  diffuse  through  the  polymer,  while  k2  may  be  different  from  k3  due  to 
blocking,  free  volume,  and  electrostatic  considerations.  In  order  to  isolate  these  two 
effects,  we  have  also  performed  simulations  with  P12  =  k i/k^  and  k2=  k3,  and  then  also 
with  P23  =  ki/k3  and  k|  =  k2. 

For  the  simulation  results  presented  in  this  study,  the  system  size  is  N=  10,000. 
Smaller  simulation  sizes  of  1M=1000,  3000,  5000,  and  7000  have  also  been  used.  The 
simulation  results  with  N=I000  differ  significantly  from  the  simulations  with  larger  N. 
On  the  other  hand,  simulations  with  larger  N  agreed  quantitatively.  This  suggests  that  the 
trends  reported  in  this  study  are  not  dependent  on  the  system  size.  Additionally,  for  the 
case  of  no  cycle  formation  and  equal  B3  reactivity,  we  compared  our  weight-averaged 
molecular  weight  with  the  analytical  theory  of  Stockmayer,14  and  the  error  is 
approximately  1%.  Clearly,  if  one  is  interested  in  describing  the  approach  to  gelation 
with  no  bound  on  the  molecular  weight,  then  the  system  size  would  also  need  to  approach 
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a  macroscopic  number  of  monomers  ( 1 and  studies  have  been  performed  to  quantify 
this  tradeoff.36  However,  in  these  simulations,  only  molecular  weights  up  to  500.000 
g/mol  are  presented,  since  this  is  the  range  of  the  experimental  data,  furthermore,  the 
error  in  the  SEC  measurements  is  approximately  10%,  so  an  error  of  1%  in  the  model 
predictions  is  not  significant. 

III.  Results  and  Discussion 

Monte  Carlo  simulations  have  been  carried  out  In  order  to  assess  the  effects  of 
cyclization,  unequal  reactivities,  and  end-capping  reactions  on  the  polymer  structure 
development.  Molecular  weights  of  As  (PPG- 1000)  and  83  (TMT)  are  1060  and  252 
g/mol.  whereas  the  molecular  weights  of  two  types  of  end-capping  reagents.  PPG-M- 
1000  and  dodecanoi,  are  1200  g/mol  and  187  g/mol,  respectively.  The  simulation  system 
containing  10,000x10,000  A  2  and  B3  monomers  yielded  molecular  weights  in  the  same 
range  as  those  observed  in  the  experiments.  Similar  to  the  experiments,  simulations  have 
been  initialized  with  an  equal  number  of  A2  and  B3  monomers  in  the  system.  The 
simulations  results  plotted  are  averages  over  50  independent  realizations, 
a.  Effects  of  cyclization  reactions 

The  melt  polymerization  was  previously  considered  to  be  sufficiently 
concentrated  so  that  the  cycle  formation  would  be  negligible.'  We  first  investigate  the 
extent  of  cyclization  via  the  simulations.  Extent  of  cyclization  (EOC)  is  defined  as  the 
fraction  of  reactions  between  A  and  B  that  is  intramolecular,  and  is  a  quantity  that  has 
been  measured  previously  using  MALDI-ToF  for  linear  polymers  but  is  less  reliably 
measured  in  hyperbranched  polymers  due  to  the  many  possible  isomers.14  ,7 


10/26 


Appendix  B 

Submitted  to  Macro  molecules 

Figures  1(a)  and  1(b)  show  the  experimental1  and  simulated  evolution  of  the 
weight-average  molecular  weight  and  the  polydispersity  PD1  as  a  function  of  A; 
conversion,  for  the  10,000x10,000  At+Bj  system  with  different  y  values.  The  reactivity 
ratio  p  is  I  and  no  monofunctional  reagents  are  present.  For  all  y  values,  a  slow  increase 
in  Mw  and  PDI  values  was  observed  until  about  80%  A;  conversion.  Above  80%  At 
conversion,  a  sharp  increase  in  Mw  and  PDI  takes  place  in  all  the  systems,  except  for  the 
one  with  the  highest  level  of  cyclization  (y  =  1).  The  experimental  data  for  PDI  and  Mw 
are  most  consistent  with  a  value  of  y  around  l O'2.  The  ideal  limit  of  no  cycle  formation, 
modeled  by  Flory,  is  the  solid  curve  with  y  =  0. 

Mw  and  PDI  development  in  systems  with  cyclization  ratios  in  the  range  of  y  =  0 
to  10‘2  all  agree  reasonably  well  with  the  experimental  data.  Due  to  the  variability  of  the 
experimental  measurements,  the  goal  of  the  modeling  is  not  to  match  the  experiments 
exactly,  but  to  assess  the  magnitude  of  the  various  effects.  The  simulation  with  y  =  10'1 
has  a  lower  Mw  and  PDI  at  high  conversions,  compared  to  the  experimental  data.  This 
suggested  that  the  extent  of  cyclization,  which  is  the  fraction  of  reactions  that  are  cycle 
forming,  was  quite  low  during  the  experiment.  Figure  2  shows  that  the  system  with  y  = 
10'1  reaches  an  extent  of  cyclization  of  0.09  at  90%  At  conversion.  Interestingly,  even 
such  a  low  extent  of  cyclization  dramatically  suppressed  the  Mw  and  PDI  as  illustrated  in 
Figures  1(a)  and  1(b).  Figure  2  also  shows  that  with  y  =  10‘3  or  y  =  10'*,  the  extent  of 
cyclization  is  less  than  3%.  This  result  supports  the  original  hypothesis,  that  melt 
polymerization  would  suppress  the  effect  of  cycle  formation  on  molecular  weight  and 
gelation. 
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Another  important  characteristic  of  hyperbranched  polymers  is  the  fraction  of 
dendritic  units  fp,  which  is  directly  proportional  to  the  extent  of  branching  in  the  system. 
The  development  of  fp  at  different  y  levels  is  shown  in  Figure  1(c)  as  a  function  of  A; 
conversion.  As  y  is  increased  from  y  =  I  O'*  to  y  =  10'!,  an  increase  in  fD  is  observed  at  A2 
conversions  of  60%  and  above.  This  trend  was  expected,  since  cyclization  reactions 
enhance  the  number  of  dendritic  groups  via  the  formation  of  small,  fully  reacted 
polymers  with  no  free  groups.  However,  the  evolution  of  fp  is  not  consistent  with  the 
experimental  data  for  any  value  of  y.  For  all  y,  the  simulations  predict  a  higher  fraction 
of  dendritic  units, 
b.  Effects  of  unequal  reactivities 

The  disagreement  of  fp  between  the  experiments  and  simulations  suggests  that  there  is  an 
additional  effect  that  suppresses  the  amount  of  branching  during  the  experiments,  other 
than  the  effect  of  cyclization  reactions.  It  could  be  attributed  to  the  lower  reactivity  of 
free  B  groups  in  linear  units  relative  to  the  free  B  groups  in  the  terminal  units  or 
completely  unreacted  B3  monomers.  However,  this  would  also  reduce  the  molecular 
weight.  In  order  to  assess  this  trade-off  quantitatively,  additional  simulations  have  been 
performed  during  which  the  reactivity  of  free  B  groups  is  modified,  based  on  the  overall 
state  of  the  B3  monomer.  In  these  simulations  presented  in  this  section,  y  =  0  since  the 
previous  section  demonstrated  that  cyclization  did  not  play  a  major  role  in  this  system. 

Figure  3  shows  the  evolution  of  the  simulations  with  different  levels  of  reactivity 
ratio  p.  With  p  =  1,  the  simulations  are  identical  to  those  in  Figure  1,  while  larger  p 
reduces  the  amount  of  branching  and  also  the  molecular  weight  and  distribution.  The 
reduction  in  molecular  w-eight  was  expected  since  the  polymers  arc  becoming  more 
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linear,  but  the  reduction  of  molecular  weight  and  polydispersity  is  not  as  dramatic  as  it  is 
with  cyclization.  Even  for  the  extreme  value  of  p=  10,  gelation  is  delayed  but  not 
completely  suppressed.  In  Figure  3(c),  the  increase  in  p  causes  a  decrease  in  fb,  as 
expected,  but  the  shape  of  the  curves  does  not  match  the  experimental  data.  In  the  data, 
the  fraction  of  dendritic  units  rises  quite  high  at  late  conversion,  but  is  very  low  at  earlier 
conversions.  The  1:1  molar  ratio  of  Aj  and  B3  is  important  in  understanding  this 
behavior.  Due  to  the  1:1  molar  ratio,  there  is  an  excess  of  B  groups,  so  at  full  A 
conversion,  only  2/3  of  the  B  groups  have  reacted.  If  the  reactivity  of  the  third  B  group  is 
strongly  reduced,  then  there  will  be  few  dendritic  units  that  form,  but  this  is  not  observed 
in  the  experiments.  The  fact  that  the  experiments  eventually  reach  a  large  value  of  fb 
near  that  predicted  with  p  =  1  instead  suggests  that  a  B  group  in  a  linear  unit  has  a  similar 
reactivity  to  a  B  group  in  a  terminal  unit. 

The  low'  values  of  fu>  around  60%  conversion  are  more  consistent  with  a 
suppression  of  the  reaction  of  terminal  units  relative  to  free  units,  as  shown  in  Figure  4 
for  various  levels  of  P12.  Recall  that  pi2  =  k]/k;,  with  k:  =  k3.  The  high  reactivity  of  the 
free  B3  could  be  due  to  its  mobility,  as  well  as  the  fact  that  B  groups  in  polymers  may  be 
partially  blocked  by  other  portions  of  the  polymer.  This  trend  in  fu  is  much  more 
consistent  with  the  experimental  data.  A  high  value  of  pi 2  =  10  best  matches  the  fb 
measurements,  while  a  somewhat  tower  value  of  P12  near  1.5-2  agrees  best  in  the 
molecular  weight  distribution.  Our  conclusion  based  on  the  simulations  in  Figures  3  and 
4  is  that  a  suppressed  reactivity  of  the  third  B  group  in  the  linear  unit  is  not  consistent 
with  the  observed  data.  The  more  consistent  explanation  is  that  the  free  B3  monomers  are 
more  mobile  and  therefore  react  faster  than  B3  in  polymer. 
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c.  Effects  of  end-capping  reagents 

The  third  and  final  simulation  study  considers  the  addition  of  monofunctional 
reagents  to  the  A2+B3  system.  Stockmayer’s  theoretical  studies  of  highly  branched 
polymers  indicated  that  addition  of  a  monofunctional  end-capping  reagent  should  shift 
the  gel  point  to  higher  monomer  conversion  values.16' ls'  w  Thus,  delaying  the  gel  point  by 
terminating  some  of  the  B  functionalities  is  the  main  motivation  behind  using  end¬ 
capping  reagents  in  this  system.  A  molar  ratio  of  monomers  as  A2:B3:E  =  1:1:1  was  used 
in  the  experiments  to  ensure  that  residual  B  end  groups,  which  would  be  expected  in  an 
A::Bj  =  1:1  system  at  full  A2  conversion,  do  not  remain  in  the  system  at  full  conversion. 

Figure  5  provides  a  comparison  of  the  change  in  the  evolution  of  Mw  for  the  addition 
of  PPG-M-1000  (1200  g/mol)  at  the  beginning  of  the  reaction.  The  A  conversion  that  is 
plotted  also  includes  the  conversion  of  the  end-capping  reagent,  since  that  is  measured  by 
NMR.  Curves  are  shown  for  the  values  of  pi2  considered  previously  in  Figure  4.  A 
primary  observation  is  that  the  addition  of  the  end-capping  reagents  has  the  larger  effect 
when  the  reactivity  ratio  is  also  targe.  Furthermore,  this  effect  is  only  observed  when  the 
end-capping  reagents  also  have  a  higher  reactivity  than  the  A  groups  in  A2  (e  »  1).  This 
might  be  the  case  for  the  dodecanol  reagent,  if  the  end-capping  reagent  has  a  higher 
diffusivity  than  the  A2  due  to  its  lower  molar  mass. 

In  the  simulations  presented  in  Figures  5  and  6,  we  set  s  =  1000,  although  the  results 
are  similar  for  e  =  10.  When  e  » 1 ,  the  simulations  predict  that  the  end-capping  reagents 
have  a  negligible  effect  on  the  polymer  structure.  By  reacting  with  the  excess  B  groups, 
they  only  add  their  extra  mass  to  the  polymer.  In  the  opposite  limit,  when  e  »  1  and  P12 
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»  1,  the  E  groups  react  quickly  with  the  free  B3  monomers,  after  which  the  EB3  units 
begin  reacting  with  A;.  Each  B3  is  thus  bonded  to  only  two  A2  monomers,  so  the 
polymers  have  a  linear  structure. 

In  the  experiments,  it  was  observed  that  the  gel  point  was  completely  suppressed 
up  to  98%  conversion  of  each  monomer,  and  the  measured  g'  contraction  factor  from 
GPC  was  more  consistent  with  a  highly  branched  polymer  (large  fD),  Our  simulations  do 
predict  a  suppression  of  molecular  weight  w'ith  the  end-capping  reagents,  but  gelation  is 
only  delayed  and  not  completely  suppressed.  At  the  extreme  value  pi2  =  10,  a  significant 
reduction  in  fD  is  also  implied,  but  at  lower  values  of  pi2,  high  levels  of  branching  are  still 
predicted. 

In  the  presentation  of  the  experimental  results,1  we  suggested  that  ester 
interchange  of  the  A2  with  the  monofunctional  reagents  might  account  for  the  observed 
reduction  in  molecular  weight.  While  this  interchange  would  cause  a  randomization  of 
the  polymer,  at  the  high  conversion  of  98%  it  does  not  provide  a  consistent  explanation 
for  the  extreme  reduction  in  molecular  weight.  At  98%  conversion,  most  of  the  A2 
monomers  freed  up  by  ester  interchange  would  have  reacted  with  another  B3  monomer. 

Other  effects  not  included  in  the  model  could  be  causing  the  suppression  of 
gelation  observed  in  the  experiments,  such  as  the  spatial  distribution  of  the  monomers  in 
the  polymer.  These  effects  could  be  exacerbated  w'hen  the  end-capping  reagents  are 
added,  since  all  B  groups  must  eventually  react,  even  those  buried  or  blocked  in  the 
center  of  the  spherical  polymer.  Possibly,  at  high  conversion,  such  groups  are  more 
likely  to  undergo  cyclization  reactions,  which  would  suppress  the  molecular  weight. 
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An  alternative  explanation  was  also  suggested  by  our  simulations.  We  observed  that  the 
simulation  results  are  extremely  sensitive  to  the  stoichiometry  near  A2:B3:E  =1:1:1.  In 
particular,  if  there  is  a  reduction  in  the  amount  of  Bj,  then  not  all  of  the  Aj  groups  will  be 
able  to  react  with  B  groups,  and  the  molecular  weight  will  be  reduced.  This  may  be  a 
particular  issue  in  our  experiments,  since  B3  loss  may  be  facilitated  by  the  nitrogen  purge 
at  our  final  polymerization  temperature  of  1 80  °C.  Due  to  the  uncertainty  in  our  final 
stoichiometry  measurements  it  is  not  possible  to  eliminate  this  effect.  The  simulations 
suggest  that  that  our  chosen  stoichiometry  of  1:1:1  is  not  a  robust  operating  point,  due  to 
the  extreme  sensitivity  of  the  molecular  weight  on  the  stoichiometry.  Figure  6  shows  a 
comparison  of  the  simulations  with  A2:B3:E  =  1:1:1  and  with  1:0. 9:1.  At  this 
stoichiometry,  90%  A2  conversion  is  the  maximum  that  can  be  achieved,  and  gelation  is 
completely  suppressed  at  full  B3  conversion. 

4.  Conclusions 

The  formation  of  highly  branched  poly  (ether  ester)s  by  the  melt  condensation  of 
an  A:  oligomer  with  a  B3  monomer  has  been  studied  using  experiments  and  kinetic 
Monte  Carlo  simulations.  The  simulations  demonstrated  that  unequal  reactivities  can 
play  an  important  role  in  the  structure  development  of  hyperbranched  polymers,  even 
when  it  has  a  little  impact  on  the  molecular  weight.  The  results  also  indicate  that  the 
presence  of  end-capping  reagents  delays  the  gel  point.  However,  the  effect  of  end¬ 
capping  agents  also  depends  strongly  on  the  ratios  of  the  various  monomers  and  their 
reactivity  ratios.  These  results  are  motivating  our  further  study  of  the  role  of  end-capping 
reagents  in  the  A2+  B3  system.  In  summary,  the  kinetic  Monte  Carlo  simulations  provide 
a  tool  for  quantitatively  assessing  the  effects  of  simple  reaction  mechanisms  on  molecular 
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structure  evolution,  enabling  the  consideration  of  a  broader  range  of  mechanisms  than 
with  analytical  models. 
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Figure  i.  Comparison  of  simulation  and  experiment1  (*).  In  the  simulations,  the 
cyclization  ratio  y  is  varied:  y  =  0  (solid  line),  y  =  10"’  (dashed  line),  y  =  10":  (dotted  line), 
y  =  10'1  (dashed-dot  line),  (a)  Weight-average  molecular  weight  Mw  (b)  Polydispersity 
index  PD1  (c)  Fraction  of  dendritic  units  fo.  Agreement  between  experiments  and 
simulations  is  not  achieved  for  at  any  value  of  y. 


18/26 


Appendix  B 

Submitted  lo  Macromolecules 


Figure  2.  Simulation  predictions  of  extent  of  cyclization.  The  cyclization  ratio  y 
varied:  y  =  I  O'3  (dashed  line),  y  =  10‘2  (dotted  line),  y  =  I  O'1  (dashed-dot  line).  For  the  y 
0  case,  EOC  =  0. 
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Figure  3.  Comparison  of  simulation  and  experiment1  (*).  In  the  simulations,  p  =  k]/k;  = 
k^/kj.  p  =  I  (solid  line),  p  =  1.5  (dashed  line),  p  =  2  (dotted  line),  p  =  10  (dashed-dot 
line),  (a)  Weight-averaged  molecular  weight  (b)  polydispersity  index  (c)  fraction  of 
dendritic  units. 
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Figure  4.  Comparison  of  simulation  and  experiment1  (*).  In  the  simulations,  p!2  =  k|/k2 
and  k2  =  k3.  p!2  =  1  (solid  line),  pi2  =  1.5  (dashed  line),  p12  =  2  (dotted  line).  p)2  =  10 
(dashed-dot  line),  (a)  Weight-averaged  molecular  weight  (b)  polydispersity  index  (c) 
fraction  of  dendritic  units. 
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Figure  5.  Simulated  evolution  with  end-capping  reagents  added  at  the  beginning  of  the 
process,  with  A2:B.i:E  =  1:1:1.  The  molecular  weight  of  E  (PPG-M-1000)  is  Mw.cap  = 
1200  g/inol.  As  in  Figure  4,  pi2  =  kj/k^  and  k2  “  k3.  pi 2  =  1  (solid  line),  P12  =  1 .5  (dashed 
tine),  pi2  =  2  (dotted  line),  pn  =  10  (dashed-dot  line),  (a)  Weight-averaged  molecular 
weight  (b)  polydispersity  index  (c)  fraction  of  dendritic  units.  (Note:  Dendritic  units  are 
calculated  here  based  on  the  number  of  A-B  reactions.  E-B  reactions  are  not  considered 
in  the  calculation  since  they  do  not  lead  to  further  branching.) 
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Figure  6.  Simulated  evolution  with  end-capping  reagents  added  at  the  beginning  of  the 
process,  with  variation  in  stochiometry:  A2iB3:E  =  1:1:1  (solid  line),  A2:Bj:E  =  1:0.9: 1 
(dotted  line  with  markers).  The  molecular  weight  of  E  (PPG-M-1000)  is  Mw.cap  =  1200 
g/mol.  pi2  =  kj/k2  =  1.5  and  k;  =  k3.  (a)  Weight-averaged  molecular  weight  (b) 
polydispersity  index  (c)  fraction  ofdendrtc  units. 
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